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SYNOPSIS 


The present thesis outlines a combined numerical and 
theoretical study of convective heat transfer in a fluid 
saturated porous medium. 

The subject of convective transport process'es in porous 
materials has received a great deal of attention in recent 
years owing to its importance in many practical applications 
such as the development of geothermal energy technology, 
petroleum reservoirs and building thermal insulation systems. 

The analysis of fluid flow and heat transfer is based on 
the transport equations and one has to consider the appropriate 
form of governing equations in the context of flow through 
porous media. The use of Darcy's law as the momentum 
conservation equation has now become common place and it is well 
known that this law does not take into account the effects of a 
solid boundary and the inertial forces. Instead of considering 
only the potential nature of the Darcy's equation, some 
researchers use the Brinkman model which also includes the 
viscous shear term while some others use Forchhe imer-extended 
Darcy model to deal with the inertial forces. There is yet 
another group of researchers who prefer to use Nav i er-Stokes 
type equation with Darcy resistance term included. 

One of the primary objectives of the present work is to 
study buoyancy induced flows in a saturated porous medium based 
on a generalized momentum equation which contains the non-linear 



convective terms as well as the viscous shear terms in addition 


to the Darcy resistance term. Numerical techniques as well as 
the analytical method is employed to obtain the solution of 
governing non-linear coupled partial differential equations. 

The first chapter is introductory and is concerned with the 
genera; background of the work presented in this thesis. 
Description and various properties of a porous medium along with 
the conservation equations of mass, momentum and energy are 
discussed in brief. Literature survey of the published work 
related to the different problems analyzed here, is presented. 

The second chapter deals with the prediction of steady 
natural convection adjacent to a semi-infinite vertical surface 
embedded in a saturated porous medium, subject to a prescribed 
constant wall temperature or the wall temperature as a power 
function of distance from the leading edge. Finite difference 
method, based on highly implicit scheme with variable mesh size, 
is employed. It is found that the vertical velocity component 
reaches a maximum a little way out from the wall and then 
decreases. The horizontal velocity component is always directed 
towards the wall. The temperature decreases monoton i ca 1 1 y away 
t rom the wail. 

Third chapter is concerned with the study of sfeady, 
two-dimensional buoyancy induced flow in a ctangular vertical 
porous enclosure with specified temperature on all the four 

m 

boundaries. The governing equations in stream function 
vorticity form are numerically solved by upwind finite 


difference method. 


A line-by-line numerical procedure is 



X 


adopted to reduce the difference equations in the tridiagonal 
form. The stream function, vorticity and temperature fields are 
obtained for various values of the dimensionless parameters. 
Flow patterns and isotherms are plotted and results on local and 
average Nusselt numbers are reported. The flow is unicellular 
filling the entire half of the enclosure. The velocities are 
higher near the wall and lower near the line of symmetry. The 
isotherm patterns are more distorted when the convection is 
domi nant . 

In the fourth chapter, free convective heat transfer in a 
vertical porous annulus whose inner wail is heated at a constant 
temperature . outer wall is isothermal ly cooled and the top and 
bottom are insulated, is i t ^/est i gated . The governing equations 
are axisymmetric form of two-dimensional continuity equation, 
generalized momentum equation and the energy equation. The 
velocity and pressure formulation is used and finite element 
method based on Galerkin approach is applied to solve the 
governing equations. Velocity and temperature fields and the 
local Nusselt number are shown as a function of the 
dimensionless parameters of the problem. Results on average 
Kjsselt number along the heated wall are also presented. 

In the fifth chapter, the problem of convective flow along 
an infinite vertical wall of constant tempe ature is analyzed. 
The far stream temperature is also constant but differs from 
that of the wall. The free stream oscillates about a constant 


mean and the effect of natural convection on oscillatory flow 



is studied. Finite difference method with Crank Nicolson scheme 


is used to predict the velocity field. Skin friction at the 
wall is also calculated. 

A variational principle based on Onsager's linear theory of 
irreversible processes is applied to study the 
thermohydrodynami ca 1 stability of natural convection in a 
horizontal porous layer bounded by two rigid walls in the last 
chapter. The critical wave numbers and the corresponding 
Ray 1 e i gh-Darcy numbers are obtained for various values of Darcy 


number . 



CHAPTER I 


INTRODUCTION 


1.1 OVERVIEW 

The transport of heat by convection in porous media is a 
process that finds frequent application in a broad spectrum of 
disciplines ranging from chemical engineering to geophysics. 
Examples of convection through porous media may be found in 
insulation systems, petroleum and reservoir engineering and 
ceramic industry. The same processes play a vital role in the 
design procedures for putting down the fire in goaf areas of 
coal mines and also in the design of pebble bed nuclear reactors 
and nuclear waste storage devices. Since convective heat 
transfer through porous media has such an overwhelming impact on 
our life, a good knowledge of these processes is necessary for 
the concerned scientists and engineers. An attempt is made here 
to study the flow of energy carrying fluids through porous 
media. The subject matter covered is the mathematical 
description of fluid flow and heat transfer through saturated 
porous media, with emphasis on the numerical solution of the 
governing partial differential equations. The improved 

numerical methods have permitted the use of general mathematical 
models capable of handling complex problems without 
approx inmt i ng the actual flow problems by a simpler one. 

Numerous studies of convection in porous media have been 
conducted using the Darcy's law which expresses that the 
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velocity is proportional to the pressure gradient. The Darcy’s 
model neglects the effects of a solid boundary and inertial 
forces. When the porous medium is bounded by impermeable 
surface, the no-slip condition must be applied and the usual 
Darcy equation is then improper to describe the flow, at least, 
near the surface. Brinkman <1947> was the first to include the 
viscous shear term in addition to the Darcy resistance term in 
the flow equation. The viscous shear stress term takes into 
consideration the boundary effect and removes the inconsistency 
of Darcy's law with the no-slip condition. The inertia forces 
become significant in the case of porous media with high 
permeabilities and also if the fluid velocity is high. To 
incorporate the inertia effects. Forchheimer <1901) considered 
additional velocity square term in the Darcy equation which 
approximately accounts for the losses in the pressure due to 
flow separation in the pores. In many applications, the inertia 
forces are treated by adding non-linear convective terms in the 
flow equations. The boundary and inertia effects not included in 
Darcy's model, may alter the velocity fields and heat transfer 
characteristics. It is seen that these effects reduce the heat 
transfer rate for a given flow driving force. 

Current research is focussed on generalizing some 
interesting flows of viscous fluid through porous medium using 
such kinds of equations that take into account the effects of 
the viscous stress as well as the inertia forces. Some typical 
cases of convection in porous media are studied while the 
ensuing system of partial differential equations are solved by 



5 


finite difference or finite 

e 1 ement 

methods 

and 

i n one 

case 

irreversible thermodynamics 

i s 

employed. 

The 

velocity 

and 

temperature distributions along 

wi th 

the heat 

transfer 

rates 


have been obtained. 

1.2 DESCRIPTKW AND PROPERTIES OF POROUS MEDIA 

It is first of all necessary to describe some of the basic 
concepts associated with the physics of flow through porous 
media. We begin with the term "porous media". A great variety 
of natural and artificial materials are porous. Soil, sand 
filters, fibrous aggregates and sandstones are examples of 
porous media. To give a precise definition of the term porous 
media is not very easy. Bear, Zaslavsky and Irmay (1968) define 
porous medium as a portion of space occupied by heterogeneous or 
multiphase matter. It should include a solid phase distributed 
throughout the medium. The solid phase is referred to as solid 
matrix and the space which is not a part of the solid matrix is 
called void space or pore space. At least some of the pores are 
interconnected and the pores which are not interconnected are 
considered as the part of the solid matrix. While the 
microscopic behaviour of a single particle or fluid in a single 
pore obeys well-known physical laws, the macroscopic behaviour 
of the porous material is described by average properties. A 
porous medium is characterized by a variety of geometrical 
properties. Following are some of the main characteristics of 


porous media. 



4 


Porosity: 

One of the most important properties of porous media is 
porosity which is a measure of pore space and hence of the fluid 
capacity of the medium. Porosity is defined as the ratio of the 
void space to the total volume and is expressed either as a 
fraction of one or in percent. 

One may distinguish two types of porosities, namely, 
absolute and effective. Absolute porosity is the percentage of 
total void space with respect to the bulk volume regardless of 
the interconnection of the pores. A porous structure may have 
considerable absolute porosity and yet have no conductivity to 
fluid for lack of pore interconnection. Effective porosity is 
the percentage of interconnected void space with respect to the 
bulk volume. It is an indication of conductivity to fluid but 
not necessarily a measure of it. 

Permeability: 

The permeability of a porous material may be defined as its 
fluid conductivity, or ability to let fluid flow within its 
interconnected pore network. Hence it is natural to expect a 
relation between permeability of a medium and effective porosity 
but not necessarily with absolute porosity. All factors, 
namely, grain size, grain packing, grain angularity, grain size 
distribution, cementation and consolidation, which influence 
effective porosity also influence permeability. 
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Specific Surface: 

The specific surface of a porous material is defined as the 
interstitial surface area of the pores per unit bulk volume of 
porous nnaterial. 

Fluid Saturation: 

There is yet another important factor to influence the 
dynamics of fluid flow through porous media — the fluid content 
of the porous structure. Saturated flow is a term used to 
describe the flow of a fluid through a porous medium where the 
entire void space is filled with the fluid. The void space of a 
porous medium may be partially filled, the remaining void space 
being occupied by air. On the other hand, two or more 
immiscible liquids may jointly fill the void space, for example, 
petroleum reservoir rocks normally contain both petroleum 
hydrocarbons and water. The saturation of a porous medium with 
respect to a particular fluid is defined as the fraction of the 
void volume of the medium filled by the fluid in question. 

All these macroscopic properties of porous media have 
significance only for samples of porous materials containing 
relatively large number of pores. A complete description of 
these bulk properties, their experimental determination 
procedures and some theories showing their interrelation are 
discussed by Muskat <1946>, Scheidegger (1957), Collins (1961) 


and Bear (1972). 
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1.3 GOVERNING EQUATICWS F(Mi CONVECTION THROUGH POR<XIS MEDIA 

The problem of complete investigation of fluid flow and 
heat transfer characteristics is to study the velocity 
distribution, the temperature distribution, the pressure and 
density of the fluid. These are functions of spatial 
coordinates and time. The path of heat and fluid flow through a 
porous medium is complex and the flow geometry changes 
continuously from one region of the material to another; for 
example, the velocity is large in small pores and small in the 
larger ones. Although the flow of fluid and energy through 
porous media is complicated, it usually is sufficient to 
determine average velocities, average temperature distribution, 
average flow paths, or the pressure distribution in the 
material. Here, the conservation equations for porous media 
flow and heat transfer are discussed in brief. 

X. Mass Conservation: 

The mass balance equation results from the consideration 
that no fluid can be produced or destroyed in any infinitesimal 
volume of the material filled with fluid. Mathematically, for a 
three-dimensional average flow the mass conservation statement 
reads , 

^ + p div q = 0 , <1.1> 

where, q is the volume averaged velocity vector, p is the 
density of the fluid and t is the time. The symbol ^ denotes 
here the substantive derivative which consists of the local 
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contribution (in non-steady flow) and the convective 

variation (due to translation), q.grad p. For an i nconpress i bl e 
fluid with p - constant, the mass conservation equation <1.1> 
assumes the simplified form, 

divq=0. (1.2) 

II. Momentum Conservation 

The momentum equation in the fluid mechanics of porous 
media originated in the nineteenth century in connection with 
the engineering of public fountains in the city of Dijon, 
France. Most of the information regarding the experiment 
performed by the French engineer Henry P.G. Darcy (1856) and his 
law in its simplest form with discussion on its limitations and 
some of its modifications may be found in the books by 
Muskat <1946). Scheidegger (1957). Collins (1961) and 

Bear (1972). Based on his experiments Darcy (1856) discovered 
that the rate of fluid flow Q through a filter bed is directly 
proportional to the area A of the bed and to the difference 
(hi - h^) between the fluid heads at the inlet and outlet faces 
of the bed and inversely proportional to the thickness L of the 
bed, to establish the Darcy's formula. 

constant x A<h.-hj>) 

Q = - . (1.5) 

L 

This form of Darcy’s law has only a very restricted use which 
can be summarized in a more useful form to state that the 
area-averaged velocity through a column of porous material is 
directly proportional to the pressure gradient established along 
the column. Denoting the area-averaged velocity by u and the 
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pressure gradient in the channel by the Darcy’s law may be 

expressed in the differential form as . 

dp 

u = constant x . <1.4) 

dx 

Subsequent research workers established that the area-averaged 
velocity of the fluid seeping through a porous material is 
inversely proportional to the viscosity of the fluid. Therefore, 
the Darcy's law for one-dimensional flow can be written as. 


K , dp , 

u = — < - -f- ) , 

p dx 


<1 .5) 


where, p is the coefficient of viscosity. The constant of 
proportionality K is called the permeability and as stated 
earlier, is a property of the porous medium. 

It is now possible to define the permeability of a porous 
medium as the volume of a fluid of unit viscosity passing 
through a unit cross-section of the medium in unit time under 
the action of a unit pressure gradient. It is determined only 
by the structure of the medium and is entirely independent of 


the nature of the fluid. Its dimensions are those of an area. 


In practical problems, the flow is rarely, rectilinear and 
neither the direction of the flow nor the magnitude of the 
pressure gradient is known. In vectorial notation, the three 
dimensional generalization of equation (1.5) including the body 
force term is. 


q 


^ <7p - p g> , 


< 1 . 6 ) 


where, q is the volume rate of flow through a unit 
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cross-sectional area of the porous medium under the action of a 
pressure gradient Vp, and g is the acceleration due to gravity. 
The quantities p and p are averaged over a region available to 
flow that is large with respect to the pore size. The 
permeability K may vary from point to point and may not be the 
same in all directions. However, it will generally suffice to 
consider the medium to be isotropic in which K is independent of 
direction. Thus the above three-dimensional generalization of 
Darcy's law is applicable to the flow of an incompressible fluid 
through a homogeneous and isotropic porous medium. Carman <19?6> 
gives an excellent exposition of the principles of flow of gases 
through porous media. The case of compressible fluid flow is 
beyond the scope of the present investigation. 

Equation <1.6) is essentially a balance of bulk viscous 
force, gravitational force and the pressure gradient and it was 
soon noticed that Darcy's law possibly can be valid only in a 
certain seepage velocity domain, outside which more general flow 
equations must be used to describe the flow correctly. In many 
practical systems, the porous medium has a high permeability and 
consequently the speed in the porous bed is not small, making 
Darcy's law inapplicable because it neglects the inertial 
forces. Forchheimer (1901) modified the Darcy's law by including 
a second order term in the velocity and proposed the relations 
of the form. 

^ + au + bu^ = 0 . (1.7) 


and later added a third order term to give. 
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dp 2 3 

^ + au + bu + cu^ = 0 . <1.8) 

dx 

Here, u is as usual the seepage velocity, p is the pressure 

(gravity is neglected), and a.b.c are empirical constants. The 
quadratic and cubic terms were added to make the equation fit 
experimental data better. When the flow is two or three 

dimensional the obvious formal generalization of Forchheimer 
extended model (Choudhary et al . 1976) may be, 

'jp - P g + g q t e| |q| q = 0 . <1.9) 

where, jq| is the absolute velocity and ^ and ^ take the place 
of a and b respectively, in equation <1,7). The Forchheimer 
equations were postulated to account for the flow separation as 
it occurs in flow through the pores over the rigid matrix. 

Nevertheless, the Forchheimer relations are valid only if the 

fluid velocity is very high. 

Deviations from Darcy's law have been observed not only at 
high flow rates but also near a solid boundary. Mathematically, 
since the order of Darcy's law is lower than that of the 
Nav i er-Stokes equation, the no-slip boundary condition is not 
required to be imposed on the impermeable boundary. This 
motivated Br i nkman< 1 947 ) to modify Darcy's law by adding the 
classical frictional terms resulting in the famous Brinkman 
model which in vectorial notation is, 

Vp - p g + ^ q - = 0 

2 

where , V is the Laplacian. 


< 1 . 10 ) 
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Attempts have also been made to study the significance of 
nonlinear convective terms in the equation of motion (Hamel 
1934, Yih 1969, Raats 1972, Yamamoto and Yoshida 1974). 
Analytical work of Yamamoto and Ivsiamura (1976) has considered 
the convective acceleration terms and frictional terms in a 
generalized equation of motion to include inertia and boundary 
effects in the Darcy's law. Following Yamamoto and 
Iwamura (1976), the momentum conservation equation for unsteady 
flow of an incompressible fluid through a homogeneous, isotropic 
and saturated porous medium may be taken to be, 

p ^ = p g - Vp + g q . (1.11) 

Equation (1.11) is often referred to as the generalized momentum 
equation for a porous medium. It reduces to the Brinkman model 
when the inertial term on the left-hand side is omitted and can 
be further reduced to ordinary Darcy's law when the term of 
viscous stress on the right-hand side is also negligible. 
Moreover, for K oo, i.e. no solid matrix present, this 
generalized momentum equation for flow through porous media is 
converted to the Nav i er-Stokes equation which is the momentum 
conservation equation for a pure fluid region. In addition to 
giving the appropriate limiting solutions as discussed above it 
is also expected that the generalized momentum equation will 
give good results in the case of highly porous media. 
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III. ENERGY CONSERVATION 

The theory of heat transfer in porous media is based on the 
existence of local thermal equilibrium between the fluid phase 
and the solid matrix. This equilibrium is required in order to 
apply local volume averaging but is not always guaranteed (Chan 
and Banerjee 1981 >. However, most of the analytical work on 
fluid flow and heat transfer through porous media has been 
carried out by assuming such equilibrium. Neglecting the 
viscous dissipation, the energy equation obtained by 
Wooding (1957) to study cellular convection in a geothermal 
reservoir is, 

pCq.VT = V, (kVT) , (1.12) 

where, T is the temperature of both the fluid and solid phases, 
C is the specific heat of the fluid and k is the effective 
therrrtal conductivity of the porous medium. For an isotropic 

the energy conservation equation 
is of the form (Cheng 1978, 

(1 .15) 

Here, a = is the thermal diffusivity of the porous medium 

P C 

and S is the heat capacity ratio of the fluid-filled porous 
medium to that of the fluid . 

The rrass , momentum and energy conservation equations 
discussed above constitute the fundamental equations for the 
study of fluid flow and heat transfer through porous media. 


homogeneous porous medium, 
under unsteady conditions 
Bejan 1984), 

S |~ + q.VT = ot V^T . 
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Many analytical investigations have been made to establish these 
conservation equations and are summarized by Cheng (1978) and 
Bejan (1964) among others. Vafai and Tien (1981) have developed 
the governing equations, alongwith an indication of physical 
limitations and assumptions made in the derivation. 

For theoretical studies of convective heat transfer through 
a porous medium, the governing equations are obtained with the 
help of Boussinesq approximation in which the density p is 
treated as a constant in all terms in the equations of motion 
except the one in the body force term. The dependence of p on 
temperature is given by, 

p^ - p = p/? <T - T^) . (1.14) 

SL gL 

where, suffix 'a' stands for the value of the quantities in the 
ambient medium and ft is the coefficient of thermal expansion. 

In natural convection flows, the body force term can be 

expressed as buoyancy term. The hydrostatic pressure with the 

body force acting on the fluid constitutes the driving mechanism 

for the natural convective flows. The local pressure p may be 

broken down to two terms, one due to the motion of the fluid, P, 

known as the viscous pressure and the other due to the 

hydrostatic pressure p . that is. 

a 

p=P+p. (1.19) 

The vector (p g - Vp) in the momentum equation (1.11) can be 
expressed as, 

pg-Vp=pg - <'7p + VP) = p g - (p_g + '^Pl . 

a a 
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where, we have used Vp = p g in the gravitational field. Using 

8l 3. 

<1.14>. finally we get. 

p g- 7p = (p - p ) g - VP = ~pft <T - T > g - VP . 

3 3 

Hence the equetion of motion (1.11) for natursl convection in a 
saturated porous medium takes the form, 

^ - g (T - T^) - 1 VP + vV^q - ^ q , <1.16> 

where, v represents the kinematic viscosity of the fluid and the 
term - g/3 (T - T ) is known as the buoyancy term. It is to be 

OL 

noted that in natural convection, the velocity and temperature 
fields are coupled. 

1.4 THE BOUNDARY LAYER EOUATIONS IN POROIS MEDIA 

The basic concepts in employing the boundary layer 
approximation to fluid flow and heat transfer through porous 
media are very similar to those in the classical viscous theory 
originated by Ludwig Prandtl in 1904 for forced flow. If a 
solid body of temperature is placed in a uniform fluid stream 
U and temperature T , then the velocity changes from zero to U 

CD ■ 00 00 

and the temperature varies from T to T in a region situated 

O 00 

relatively close to the solid body. The thin layer in which the 
velocity increases from zero to its full value is called viscous 
boundary layer and in a similar way, the temperature varies from 
T to T in a thermal boundary layer. Boundary layer analysis 

O 00 

is applicable to natural convection flows also in which the 
pressure in the region beyond the boundary layer is hydrostatic 
unlike in the case of forced flow where it is being imposed by 
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an external flow. According to the boundary layer concept, 
therefore, the flow field is divided into two regions, first the 
slender boundary layer in which convection of momentum and 
energy occur and the other, outside the boundary layer where 
inviscid flow prevails. Such a division of the flow field 
brings about a considerable simplification to the mathematical 
theory of the motion of fluid of low viscosity. The boundary 
layer equations are derived by comparing the magnitude of the 
various terms in the momentum and energy equations. Certain 
expressions in the differential equations governing the flow are 
neglected due to their relative smallness as compared to those 
reta ined. 

The first paper dealing with the boundary layers in porous 
media appears to be that of Wooding (196?) to treat the problem 
of free convection about a line source and a point source, as 
well as for free convection about two finite heated vertical 
plates embedded in a porous medium. Wooding pointed out that 
when the dimensions of a convective system in a saturated porous 


medium are sufficiently great 

, diffusion 

effects can 

be 

neglected 

except in regions 

where the 

gradients of 

fluid 

propert i es 

are very large. 

A boundary 

layer theory 

was 

deve 1 oped 

for vertical plane 

flows in such regions. 

The 


boundary layer approximations were also invoked by McNabb <1965) 
who studied free convection flows in a saturated porous medium 
above a horizontal heated plate. The book by Bejan <19841 
contains an excellent coverage on boundary layers in porous 


media. 
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Making the usual boundary layer assumptions, the 
generalized momentum equation for a steady two-dimensional 


incompressible fluid flow takes the form. 


2 

, du . dp . 3 u 

p (u -3— + V -5—) = - + p , 

ffx 0y dx - 2 

3y 

and the energy equation becomes. 
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where, u and v are the velocity components in x and y directions 
and T is the temperature. For natural convection the 
two-dimensional momentum equation in steady state assumes the 
followirg form. 
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1.5 BACKGROUND OF THE PRESENT WORK 


A considerable amount of work has been published in the 
past five decades on thermal convection in porous medium. Early 
research credit goes to Horton and Rogers <194/'>, Morrison et 
al. <1948), Rogers and Morrison (1950) and Rogers et al . <1951). 
These authors considered this kind of flow in attempting to find 
the distribution of salt in subterranean sand layers. The 
problem of setting up of convection currents to discuss the 
possibility of convective flow in a porous medium was 
investigated theoretically by Lapwood (1948) and confirmed 
experimentally by Katto and Masuoka (1967). Eversince, this 
field has never lost its momentum because of its wide 
applications in science and technology. Wooding <19^7-1963) and 
Elder (1967) made significant contributions to this area in 
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connection with the studies of transport processes occurring in 
geothermal reservoirs. The basic development of the use of 
boundary layer approx i rra.t i on to examine the convective flow 
through porous media was done by Wooding (1965) and 
McNabb (1965). Much of the previous work on heat transfer in 
porous media and geothermal systems is to be found in extensive 
review article by Cheng (1978). 

Numerous investigations have been devoted to the natural 
convection in a porous medium adjacent to flat plates. The 
first investigation on free convection flow about a vertical 
flat plate embedded in a porous medium using boundary layer 
approximation was reported by Cheng and Minkowycz (1977). who 
obtained similarity solutions with the assumption that wall 
temperature is a power function of distance from the origin. 
Cheng<1977) revealed that solution to problems of prescribed 
heat flux can be obtained by a simple transformation of 
variables. Using an order of magnitude estimate Cheng (1978) 
obtained the boundary layer equations for a porous layer 
adjacent to a heated vertical surface and analyzed the problem 
using integral method with different assumed profiles. The 
local Nusselt number obtained from the integral method appeared 
to be in good agreement with that of the similarity solution. 
Cheng and Chang (1979) applied perturbation method to the 
problem of natural convection adjacent to vertical or horizontal 
heated surfaces. They assumed the power law variation in the 
wall temperature and solved perturbation equations of first 
order to give boundary layer solutions. The numerical solution 
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of the boundary layer equations was given by Na and Pop <198>> 
to predict the steady natural convection heat transfer from an 
impermeable vertical surface in a saturated porous medium 
subject to a prescribed non-uniform wall temperature or to a 
prescribed non-uniform wall heat flux. Higher-order effects of 
Darcian free convection boundary layer flow adjacent to a 
semi— i nf i n i te vertical flat plate with a power law variation of 
wall temperature were examined theoretically by Cheng and 
Hsu <1984) using the method of matched asymptotic expansions. 
Natural convection from vertical plates in porous media was 
analyzed by Ingham and Pop <1987) for cases where the leading 
edge of heating plate is at an arbitrary distance above an 
impermeable horizontal boundary and the plate is maintained at a 
constant temperature or uniform heat flux. Gov i ndara j ulu and 
Malarvizhi <1987) considered the problem for various wall 
temperatures and injection velocity conditions to obtain 
approximate series solution. 

Boundary layer analysis for natural convection in a fluid 
saturated porous medium, adjacent to horizontal surfaces with 
wall temperature being a power function of distance from the 
leading edge, was performed by Cheng and Chang <1976>. 
Similarity solutions for the convective flow above a heated 
horizontal impermeable surface or below a cooled horizontal 
impermeable surface are obtained and applications to convective 
flow above the heated bedrock or below the cooled caprock in a 
liquid-dominated geothermal reservoir were discussed. Free 
convection about an impermeable horizontal surface in a porous 
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medium with prescribed wall temperature was considered by Chang 
and Cheng (1985) and solutions were obtained by asymptotic 
matching procedures. Chandrasekhara et al . <1984) also found 
similarity solutions for free convection about an impermeable 
horizontal surface in a saturated porous medium of variable 
permeability. The results of the analysis had a bearing on 
convective flows in geothermal reservoirs. Natural convection 
for horizontal surface in porous media bounded by an impermeable 
vertical wall was examined by Ingham and Pop <1987) who applied 
the method of matched asymptotic expansions for n«.themat i cal 
treatment . 

Little attention has been given to transient convection in 
a porous medium. Johnson and Cheng <1978) carried out a 
systematic research of similarity solutions for free convection 
about flat plates in porous media for steady and transient 
situations. Similarity solutions were shown to exist for steady 
free convection about vertical surfaces when the temperature 
difference between the wall and the environment varies either in 
power-law or exponential forms and for horizontal surfaces 
according to power law form. Several very specific solutions 
for unsteady free convection about flat plates in a porous 
medium were shown to exist. The problems of transient free 
convection in a porous medium adjacent to a vertical 
semi - i nf i n i te fiat plate with a step increase in wall 
temperature and surface heat flux were considered by Cheng and 
Pop <1984). Ingham et al. <1982) presented solutions of the 
boundary layer equations which describe the transient free 
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convection flovg in a saturated porous medium near a vertical 
flat impermeable surface which is suddenly cooled to the ambient 
fluid temperature. 

Much work has been done on the study of convection in 
enclosed porous medium bounded by solid plane surfaces. 
Theoretical and experimental studies of natural convection in 
confined porous media were carried out by Holst and Aziz (1972). 
Critical Rayleigh number for convection in a rectangular box 
filled with fluid saturated porous material heated from below 
was determined by Beck (1972). Natural convection in a vertical 
porous layer of large temperature difference between the 
vertical walls was studied theoretically by Weber (1975) with 
the assumption that boundary layers develop along the vertical 
walls. Bejan (1979) modified the Weber’s theory for convection 
in a vertical porous layer by combining Weber’s boundary layer 
solution with zero energy flux boundary conditions at the top 
and bottom walls. Convection in a cavity enclosing a Darcy 


medium 

driven by 

heating in 

the 

hor i zonta 1 

was analyzed 

by 

Wa 1 k e r 

and Homsy 

(1978) where 

the 

solut ions 

are governed 

by 


dimensionless parameters of Darcy-Rayleigh number and the cavity 
aspect ratio. Bejan and Tien (1978) have reported an analytical 
method to study the corner effects in shallow cavities. Steady 
convective motion inside a rectangular cavity filled with a 
porous material were examined by Simpkins and Blythe (1980) in 
the large Rayleigh number limit for flows driven by horizontal 
temperature gradient. An integral relations approach was used 
to determine the boundary layer structure on the side walls. 
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Numerical solutions for horizontal porous layer through the 
Galerkin form of the finite element method are due to Hickox and 
Gartling <1981). Heat transfer rates were presented for aspect 
ratio ranging from 0.1 to 0.9 and Rayleigh numbers in the range 
29 to 200. Heat transfer results for low Rayleigh number are 
also due to Bankvall <1974) and Burns et al . (1976), both via 
numerical solutions. 

The results of a numerical simulation of natural convection 
in a vertical rectangular porous enclosure subjected to a 
horizontal temperature differential were presented by Shiralkar 
et al . (1985) for substantially larger values of the 
Darcy-Rayleigh number. The effects of aspect ratio on flow 
behaviour and heat transfer rates in a rectangular porous cavity 
have been reported by Prasad and Kulacki <1984a). In another 
paper Prasad and Kulacki <1984b) presented numerical results for 
two-dimensional, steady. free convection for a rectangular 
cavity with constant heat flux on one vertical wall, the other 
vertical wall being isothermally cooled. The numerical solution 
of Prasad (1987) deals with convection in a rectangular cavity 
filled with a heat generating Darcy porous medium. 

Recently natural convection in vertical annuli has also 
been studied in the presence of rigid matrix in Darcy’s regime. 
By using a finite element numerical method. Hickox and 
Gartling (1982) have obtained heat transfer results for tall 
annulus with the inner wall heated at a constant temperature and 
the outer wall isothermally cooled, the top and bottom being 
insulated. Similar problem was considered by Havstad and 
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Burns (1982) who used finite difference method, perturbation 
technique and an approximate analysis for investigation. The 
results obtained by the above authors are applicable for only 
low Rayleigh number heat transfer. The work of Philip <1982) 
also describes free convection at small Rayleigh number. Prasad 
and Kulacki (1984c) considered the same thermal boundary 
conditions but for larger ranges of Rayleigh number and radius 
ratio. A numerical and experimental investigation of natural 
convection in saturated porous media in the annulus between 
vertical concentric cylinders with small aspect ratio was 
reported by Prasad and Kulacki (1989). As the radius ratio 
(outs ide/ i ns i de) was increased, smaller aspect ratios were 
required to introduce multicellular flow structures. 
Prasad (1986) numerically examined the convective heat transfer 
in a vertical annulus when its inner wail is heated by applying 
a constant heat flux. 

In describing heat transfer in porous media, many studies 
have dealt with non-Darcy flow model. The buoyancy induced 
boundary layer adjacent to vertical wall was analyzed by Plumb 
and Huenefeld (1981) using a velocity square term in the 
momentum equation. Deviations from Darcy's law were shown to 
occur for larger values of modified Grashof number. Forchheimer 
model was also used by Bejan and Poulikakos (1984) to study the 
intermediate and non-Darcy flow regimes for natural convection 
along vertical plate embedded in a porous medium. Lai and 
Kulacki (1987) studied steady non-Darcy convection from a heated 


horizontal surface. 


Inertial (non-Darcy) effects at high 
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Rayleigh number in a vertical porous layer subjected to uniform 
heat fluxes at the sides were numerically evaluated using the 
Forchheimer model <Poulikakos 198?) . Numerical experiments in a 
porous layer heated from the side were carried out by Poulikakos 
and Bejan <1983> to confirm the non-Darcy flow prediction of an 
analytical solution and to illustrate the regime between Darcy 
flow and flow affected by inertia. Nield and Joseph (1989) 
demonstrated that the effects of inertia are significant in 
natural convection at realistically large values of the Rayleigh 
number in thin layers of saturated porous media with small 
Prandtl number. Inertia effects on buoyancy-driven flow and 
heat transfer in a vertical porous cavity were examined by 
Prasad and Tuntomo (1987). A Forchhe imer-extended Darcy model 
was shown to give superior predictions in the inertial flow 
regime . 

Several works have been published by using the Brinkman 
extended Darcy model where the viscous dissipation term is 
included in the momentum equation. The Benard problem from 
viscous flow limit to Darcy's limit has been considered by Katto 
and Masuoka (1967), as well as by Malker and Homsy (1977). These 
investigators used Brinkman model as the momentum equation. 
Rudraiah et al . (1980) and Rudraiah and Balachandra Rao (1985) 
also used the same model to understand the onset of convection 
in porous medium made up of sparse distribution of particles. 
Using matched asymptotic solutions for a three-layer model for 
natural convection at a semi-infinite vertical plate in 
saturated porous media , the Brinkman model was incorporated 
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(Hsu and Cheng 1985). The boundary effect was shown to be 
important at larger Rayleigh numbers and higher Darcy numbers. 
The Brinkman model was employed by Chan et al . (1970) for 
analyzing heat transfer in a rectangular porous enclosure filled 
with gas. Tong and Subramanian (1985) have presented a boundary 
layer analysis for the Brinkman model to investigate the effect 
of the no-slip boundary condition in rectangular enclosures 
containing a porous medium. They reported significant 
contributions of the diffusion term at high Rayleigh and Darcy 
numbers. Tong and Orangi (1986) and Sen (1987) studied natural 
convection in porous cavities incorporating the viscous terms in 
the equation of motion. 

Although as early as 1976. Yamamoto and Iwamura considered 
the transport and viscous terms in a generalized equation of 
motion, only a limited effort has been made to study the 
significance of these terms for convection in porous media. 
Two-dimensional free convection flow through a porous medium 
bounded by a vertical infinite surface was analyzed by using the 
generalized momentum equation for steady state case (Raptis et 
al . 1981) and also for unsteady state situation (Raptis and 
Perdikis 1985). Raptis (1983) and Singh et al . (1986) used this 
model to study unsteady free convective flow through a porous 
medium bounded by an infinite vertical plate, when the plate 
temperature is oscillating with time about a constant value. 
Raptis and Perdikis (1985) examined the effects of free stream 
oscillations on this geometry and Raptis and Singh (1985) 
considered the case of moving vertical wall. Very recently. 
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Lauriat and Prasad <1987) made a dimensional analysis of the 
Brinkman-extended Darcy formulation, with the transport term in 
the equation of motion to study the effects of transport and 
viscous terms on natural convection in vertical porous cavity. 

It should be noted that attempts to consider the transport 
and viscous terms for analyzing convection in porous medium are 
limited. Recent investigations suggest the incorporation of 
these terms in the momentum equation. The work in this thesis 
is devoted to the study of some convection problems in porous 
media using the generalized momentum equation. which contains 
nonlinear transport term and a term of viscous stress also. We 
consider the equation with and without boundary layer 
approximation depending on the nature of the problem. The 
results are generated by solving the complete two-dimensional 
coupled governing equations by finite difference or finite 
element numerical methods. Irreversible thermodynamics is 
employed to the classical Benard-Darcy problem. 

The problem of natural convection heat transfer from an 
impermeable semi - i nf i n i te vertical surface embedded in . a 
saturated porous medium subject to a prescribed constant wall 
temperature or the wall temperature as a function of position 
from the leading edge is considered first. The analysis is 
based on boundary layer approximation. A numerical marching 
technique based on highly implicit finite difference method with 
variable mesh size is applied for the solution. The velocity 
and temperature fields are obtained together with the rate of 


heat transfer. 
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The next problem chosen deals with natural convection in an 
enclosed porous medium bounded by solid plane surfaces. The 
geometry is limited to a rectangular one with no-slip and 
isothermal boundary conditions on all four surfaces. The 
stream-function and vorticity formulation is used and an upwind 
finite difference scheme is employed for discretization of the 
vorticity and energy equations. The results are presented in 
terms of theoretical streamlines and isotherms for the 
dimensionless parameters, the Ray 1 e i gh-Darcy number, the Darcy 
number, aspect ratio and the top wail temperature. Results are 
compared with those obtained by finite element method based on 
Galerkin approach. 

In the fourth chapter, free convective heat transfer in a 
vertical porous annulus whose inner wall is heated and the outer 
wall is i sothermal ly cooled, the top and bottom being insulated, 
is investigated. The governing equations are axisymmetric form 
of two-dimensional continuity equation, generalized momentum 
equations and the energy equation. The velocity pressure 
formulation in which the generalized momentum equations are 
considered as they are, is chosen. The velocity components, 
pressure and temperature are all obtained simultaneously. at 
each node by finite element method. The algorithm developed 
here is very general and is suitable for a wide range of 
problems . 

In the fifth chapter, transient convective flow through a 
porous medium is considered. The buoyancy effect on fluctuating 
flow through a porous medium bounded by an infinite vertical 
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wall is studied. The free stream velocity oscillates about a 
non-zero mean and the surface absorbs the fluid with a constant 
velocity. The temperature of the wall is constant and it 
differs from free stream temperature. The governing equations 
are solved by Crank-N i col son finite difference method and the 
results are obtained for various values of Grashof number, the 
permeability parameter and the amplitude parameter. 

In the last chapter of the thesis, a variational principle 
based on Onsager's linear theory of irreversible processes is 
used to study the thermohydrodynami ca 1 stability of natural 
convection in a porous medium bounded by two horizontal walls. 
The balance equations of mass, momentum and energy are 
considered to study the natural convection in a thin horizontal 
layer of saturated porous medium heated from below. The problem 
is treated by linear perturbation analysis when both the 
bounding surfaces are rigid. The critical wave numbers and the 
corresponding Rayleigh-Darcy numbers are obtained for various 
values of Darcy number. 
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CHAPTER II 


NATURAL CONVECTION ABOUT A SEMI INFINITE VERTICAL FLAT PLATE 

IN A PCM^OUS MEDIUM 


2.1 I NTRODUCTI ON 

The study of natural convection adjacent to heated vertical 
surfaces embedded in a fluid saturated porous medium has been 
the subject of several recent investigations- The introduction 
of boundary layer approx inat i on in porous medium literature 
intensified the research efforts on analytical studies of 
convective heat transfer about heated surfaces. Within the 
framework of boundary layer approximation similar to those 
invoked by Wooding (1963) and McNabb (1963), Cheng and 
Minkowycz (1977) obtained the analytical expression for boundary 
layer thickness and heat transfer coefficients to study 
convective heat transfer about an isothermal dike intruded in an 
aquifer. Johnson and Cheng (1978) analyzed free convection 
boundary layers in porous medium adjacent to flat plates 
including unsteady cases and thermal stratification. A 
comprehensive review of the early literature has been presented 
by Cheng (1978). The earlier works rely heavily on Darcy's law. 
Plumb and Huenefeld (1981) used the Forchheimer extended Darcy 
model and Hsu and Cheng (1983) used Brinkman model to study the 
problem. With few exceptions, most of the earlier theoretical 
studies are devoted to the problem of finding similarity 


solutions for steady 


free convection about the vertical 
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surfaces. In all these analysis power form variations of 
wall temperature distribution are assumed. A numerical solution 
of the problem has been obtained by Na and Pop <1985). 

The objective of the present analysis is to predict the 
steady natural convection heat transfer from an impermeable 
vertical surface in a saturated porous medium subject to a 
prescribed constant wall temperature or the wall temperature as 
a function of position from the leading edge. The generalized 
momentum equation in which the convection and viscous terms are 
taken into account has been considered. Several studies have 
been made on the basis of this generalized momentum equation 
(Singh et al . 1986. Raptis et al , 1987 and Raptis and Takhar 
1987 ). The method of solution, applied here, is close to that 
used by Hornbeck <1961), who has developed a numerical marching 
technique based on highly implicit finite difference method with 
variable mesh size. Thus, the numerical solution of the 
boundary layer equations are obtained to the required accuracy. 

2.2 BASIC ECft/ATIC»fS AND BOUNDARY CONDITICWS FOR CONSTANT WALL 

TEMPERATURE 

A two-dimensional steady flow through a fluid saturated 
porous medium occupying a semi - i nf i n i te region of the space 
bounded by a vertical impermeable surface is considered. The 
X-axis is taken along the surface and the Y-axis is normal to 
the surface. The boundary layer and Boussinesq approximations 
are assumed and the properties of the fluid and the isotropic 
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porous matrix are constant. The governing equations with 
inertia and viscous terms are. 
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where, U and 

V 

Rre the 

velocity 

components 

in the X and 


Y-directions respectively, T is the temperature, v is the 
kinematic viscosity, g is the acceleration due to gravity, ft is 
the coefficient of volume expansion of the fluid, a is the 
thermal diffusivity and K is the permeability of the porous 
medium. The subscript ao denotes the condition near the edge of 
the boundary layer. The boundary conditions are. 


X = 0; U=0. T= 

Y = 0: U=0. V= 

Y -► 00: U->0. T-> 

The equations <2.1> - <2.5> 
are put in dimensionless 
trans format ions. 
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where. L is the characteristic length in the X-direction, Da is 
the Darcy number, Gr is the Grashof number and Pr is the Prandtl 
number. In terms of the new variables, basic equations become. 
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cond i t i ons 
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2.3 SOLUTIONS BY FINITE DIFFERENCE METHOD 

Equations (2.6) - (2.8) are now written in finite 
difference form by superposing a two-dimensional rectangular 
mesh on the flow field. Step sizes Ax and Ay are taken 
respectively in the x and y directions. Subscripts <i,j> 
indicate positions in the (x. y) directions respectively so that 
i = 0 corresponds to x = 0 and j = 0 to y = 0 ; a positive change 
Ax in the x coordinate increases i by 1 and a positive change 
Ay increases j by 1 . Let the plate be represented by j = 0 and 
the edge of the boundary layer by j = n+1 . The finite 
difference equations are written as (Hornbeck 1961), 
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The difference form stated here is highly implicit, that is, not 
only all y-der i vat i ves are evaluated at <i+l) but the 
coefficients of non-linear convective terms are also evaluated 
at <i+l). The discretization is of second order accuracy in Ay 
for u and & and is suitable for the marching procedure along the 
x-direction in conformity with the parabolic nature of equations 
<2. 7) and (2.8). An iterative scheme is required to solve the 
non-linear difference equations. The difference equations 
(2.10) - (2.12) are rewritten in the linearized form as. 


k+1 _ k+1 k+1 

’^i+l.j+1 ,^i.j + l + ? + l + = 0 

Ax ' Ay 


( 2 . 1 ?) 


k+1 _ k+1 k+1 

k '^i + l,i ^i.j . ^k ^i+l,j + l ~ ^i+1,3-1 

'^i+1,3 Ax i+1. 3 2(Ay) 


k+1 j, k+1 k+1 

^i+1 , j-1 “ ^^i+1. i ^i+1, 3+1 ^ ^k 

<Ay)2 ■ ' 


u 


k+1 
i + 1 . 3 
Da 


(2.14) 



40 


gk+l _ k+1 _ k+1 

k + 1 ^i + l.j i.i , i + l.j + 1 i + 

^1+1 , j 'Ax i+1 . j 2(Ay) 


k+1 

1 i+l.j-1 


1 +1 . ) +1 • j+i 

<Ay)^ 


<2. 15) 


The superscripts k and <k+l> indicate values obtained on kth and 
<k+l)th iteration respectively. 

Equation <2.14) is. finally, written in the simplified form. 
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Similarly equations <2.1?> and (2.15) are rewritten in more 
useful forms as. 
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The linearized difference form (2.16) of the momentum equation 

written for j = 1(1 )n leads to a set of n linear equations for n 

unknown values of u at the x-location i+1. This linearized 

tridiagonal set is solved by Gaussian elimination method and 

then the solution of (2.17) for v at the x-location i+1 is 

obtained in a stepwise manner working outward from the plate. 

Once values of u and v have been determined at i+1, the 

discretized energy equation (2.18) written for j = l<l)n leads 

to a tridiagonal set of n linear equations that can be easily 

solved for the n unknown values of © at the location i+1. 

Equations (2.16) - (2.18) are solved iteratively until the 

values obtained on a given iteration agree with those obtained 

on the previous iteration within some predetermined accuracy. 

For the convergence of the iterative process underrelaxation is 

k k k 

employed. The quantities u... v.^, . and©.., . appearing in 

1+1,3 1+1.3 1+1.3 

the difference equations (2.16) - (2.18) are modified after each 
iteration as follows. 
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where, \ is 

a 

re laxat i on parameter 

and va 

lues 

of X 

i n the 

range 0 < X 

< 1 

correspond to underre laxat i on . In 

the 

present 

cal culat i on 

underrelaxation factor 

is found 

to be 

0.6. 

This 


completes the solution at the x-location i+1. A repetition of 
the above steps allows to march in upward direction along the 
wall. 

The computer program takes care of the fact that the 
velocity and temperature profiles vary rapidly near the wall and 
hence a fine mesh is required near the wall. A relatively 
coarse mesh has been used in the regions away from the wall. If 
Ay is uniform across the boundary layer thickness, the number of 
simultaneous equations increases excessively. Solution of such 
large number of equations not only requires excessive computer 
time but also involves large round-off error. Non-uniform Ay 
therefore has been used. This requires the technique of 
combining large and small mesh sizes. 

If the mesh size is to be changed from Ayj to a larger 
value Ay^ at the y location j = p as shown in Figure 1, the 
central difference form for the first or second order derivative 
at j = p requires the values of the dependent variable at a 
fictitious point j = q which is Ay^ above the point j = p. A 
dependent variable Qj+2’q point j = q is determined by 
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passing a parabola through the values Q Q 

1+1, p-1 1+1, p 

Q.+l and then interpolating for q yielding. 
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The derivatives at j = p are approximated as. 
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where, Q is given by equation (2.20). 

1 T i f q 

At the points of step size change, equations (2.21) and 
(2.22) are used in the momentum and energy equations and at 
these points equations are rewritten, as they involve central 
differences in transverse direction. This modification is not 
required in the continuity equation because only first order 
forward differences in the transverse direction are involved. 
It is ensured that the proper mesh size is used in the 
appropriate region. 

The local Nusselt number in the non-dimensional form is 
given by. 

Be , 

Nu = - ^ n • 

X By 'y=0 


(2.25) 
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This has been calculated numerically by the following difference 
formula. 


i41 .0 


i + 1 .1 
2<Ay> 


i+1 .2 


<2 . 24) 


2.4 FREE CONVECTION WITH WALL TEMPERATURE A FUNCTION OF 
PCmTION 

In section 2.2. the problem has been formulated for the 
case of constant wall temperature. The problem must be 
reformulated if the wall temperature is a function of position. 
The governing differential equations remain the same as (2.1) 
(2.3). The boundary conditions in this case are. 

X = 0: U=0. T=T . 

OD 

Y = 0 : U = 0, V = 0, T = T <X> . (2.23) 

w 

Y->cd: U-*0. T-».T . 

00 
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Introducing the dimensionless variables <2.26), the governing 
partial differential equations and the corresponding boundary 
conditions take the following form. 
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equations in the dimensionless form remain the same as <2. 6) 

(2.8) in the case of constant wall temperature. ©^(x> may be 

any function of x relevant to practical and engineering 

applications. For the sake of convenience, in the present 

investigation, 9 <x) is defined as, 

^ w 


9 <x> = x" , 

w 

where, n is a real number. The fundamental steps for solving 
the differential equations are similar to those specified in 


Section 2.3. 
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2.5 RESULTS AND DISCUSS! CW 

No quantitative statement can be made about the choice of 
step sizes. In the regions of more rapidly changing velocity 
and temperature profiles a fine mesh size is required. Since 
implicit formulation is used there is no restriction on the 
selection of step sizes Ax from the point of view of stability. 
The computations are carried out on DEC-1090 computer for Ax = 
0.02 and variable mesh size employed in y-direction are as 
f ol lows , 

Ay = 0.02 for first 50 steps. 

Ay = 0.05 for next 150 steps. 

Ay = 0.1 for next 5 steps. 

These step sizes have been finally selected by approximating the 
thickness of the boundary layer formed which corresponds to 8.0 
in the present calculations. The underrelaxation factor is 
chosen to be 0.6. Different step sizes and underre laxat i on 
factors were used and the above mentioned values were found to 
generate the best results with minimum CPU time. 

All the values of v are found to be negative which indicate 
that the horizontal velocity component is always directed 
towards the wall. The vertical velocity component u reaches a 
maximum a little way out from the wall and then decreases. The 
temperature 6 decreases monotonical ly away from the wall. 

In Figure 2 the vertical velocity field is plotted against 
y when the wail temperature is constant. The value of Prandtl 
number Pr is taken as 0.75 and that of x as 0.2, the curves are 
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plotted for different values of Da. It is observed that 
u-velocity increases significantly as the value of Da increases. 

The value of temperature at a point increases as the Darcy 
number Da decreases. However, this effect has not been shown in 
a pictorial diagram. 

The ranee of n for which the problem is physically 

consistent is considered. To this end the method given by Cheng 

and Minkowycz (1977) is followed. Since the wall temperature 

begins to deviate from T at x = 0 and convective flow must 

00 

start at this point, horizontal velocity and boundary layer 

thickness must therefore increase or at least remain constant 
with X. This implies that the exponent of x must be greater 

than or equal to zero, which results in the limit 0 ^ n S 1. 

Figure 5 shows the variation of u - velocity with y for 
different values of n. which is the exponent of wall 

temperature. The values of n are chosen to be 0.0, 0.25, 0.9 
and 0.79. The velocity at any point is found to be decreasing 
with the increase in n. 

Figure 4 shows the temperature profiles for different 

values of n. These curves are plotted for n = 0.0, 0.29, 0.9 
and 0.79 and for fixed values of x, Pr and Da. The value of 

temperature at a point decreases with n. The value n = 0.0 

corresponds to the case of constant wall temperature. 

It is customary that when a numerical method is applied to 
a problem, the results so obtained are compared with the 

available numerical solution or exact solutions to determine the 
error involved in the results. Since non-similar results are 
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not available in the existing literature for the problem 
considered in this investigation, the numerical solutions are 
obtained by setting Da oj and n = 0.0. In this case the 
governing equations correspond to free convection flow along 
vertical wall of constant temperature in the absence of porous 
matrix. The results of local heat transfer calculated for 
Pr = 0.755 are compared with those given by Schlichting (1960) 
? id theagreement is seen to be very good. It is found that the 
cal Nusselt number for x = 0.5 and 1.0 are 0.4559 and 0.5654 
by the present method while it is 0.4272 and 0.5592 respectively 
as reported by Schlichting (1968). 
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CHAPTER III 


NATURAL CONVECTION IN A RECTANGULAR POROUS 
ENCLOSURE WITH HEATED WALLS 


3.1 INTRODUCTION 

Natural convection in fluid saturated porous enclosures has 
received considerable attention for its application in a variety 
of engineering problems. A number of analytical, numerical and 
experimental studies have been performed for natural convection 
in vertical rectangular porous enclosures with isothermal 
vertical walls and adiabatic horizontal walls. The assumption 
of insulated horizontal walls may be of fundamental interest but 
there are some physical situations in which vertical heat 
transfer through horizontal boundaries, is not negligible. 
Investigations that focus on natural convection in rectangular 
enclosures in which the side walls, top wall and the bottom 
wall are at constant but different temperatures, are few in 
number. This condition permits heat transfer across all the 
boundaries, which results in considerable variation of the flow 
field. Examples of such situations include geothermal 
applications and coal mine fires in the goaf areas. 
Christopher <1986> has analyzed transient natural convection for 
the case when all the four walls are isothermal. 

Furthermore, in most of these studies it is assumed that 
Darcy’s law is applicable. Analytical studies by Weber <1975), 
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Walker and Homsy (1978). Bejan (1979) and Simpkins and 
Blythe (1980) and numerical studies by Holst and Aziz <1972), 
Hickox and Gartling (1981), Shiralkar et al. <1985), Prasad and 
Kulacki (1984) and Christopher (1986) are based on Darcy's law. 
Chan et al. (1970), Tong and Subramanian (1985) and Tong and 
Orangi (1986) used the Brinkman-extended Darcy model, which 
included the viscous force terms in their analysis and the 
no-slip boundary conditions were satisfied. Pouiikakos and 
Bejan (1985) and Prasad and Tuntomo (1987) have analyzed the 
inertia effects by using Forchhe imer-extended Darcy equation of 
motion and the significance of Forchhe imer ' s velocity square 
term are examined. Recently, Lauriat and Prasad (1987) have 
performed a numerical study for a generalized momentum equation 
and have included viscous force terms and inertial force terms 
simultaneously. All the above investigations except the one by 
Christopher (1986) describe convective heat transfer in 
rectangular porous enclosures with two adiabatic walls. 

The present investigation is concerned with the steady, two 
dimensional buoyancy- i nduced convection in a vertical porous 
rectangular enclosure with specified temperature on all the four 
surfaces. The boundaries are impermeable and the results are 
based on inclusion of boundary and inertia effects. The 
pressure field appearing in the momentum equation is eliminated 
and the stream function vorticity equation together with the 
energy equation is solved by the upwind finite difference 
method. This method leads to a simple and straight forward 
analysis and possesses a highly desirable property of numerical 
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stability. The stream function, vorticity and temperature 
fields and the local and average Nusselt numbers are obtained 
for four dimensionless parameters. namely. the Ray 1 e i gh-Darcy 
number. Darcy number, aspect ratio and the top wall temperature. 


3.2 PROBLEM STATEMENT 

The physical system shown in Figure 1 is a two dimensional 
fluid saturated porous medium enclosed by four isothermal and 
impermeable walls. The two vertical side walls are maintained 
at the same temperature T^ and the top and bottom horizontal 
walls are kept at two different temperatures T^ and T^^ 
respectively. Subject to the Boussinesq approximation, the 
continuity, momentum and energy equations for the flow of a 
Newtonian fluid are. 


^ ^ = 0 

ax ay 


<5. 1) 


,,, au ^ ., du, ap ^ , a^u ^ a^u 
ax ^ ay^ ~ ' ax ^^z 


> - 


K 


(5.2> 


/«. dv. _ 9P , . a^w , a^y . 

ax * ^ 57 ^ ~ ~ ay ^^z ^ ^^z ^ 


,, dT ^ dT , a^T , a^T 


- ^ + pg^<T-Tg> . 


( 5 . 5 > 

<5.4> 


The boundary conditions for the problem are, 
X=0:U=0, V=0. T=T , 


X = A: U = 0. V = 0. T = T, 
Y=0:U=0. V=0. T=T 


<3.5> 


b • 


= 0. V = 0. T = T 


t ’ 


Y = B: U 
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where. U and V are the horizontal and vertical velocity 
components in the X and Y directions respectively, P is the 
pressure, T the temperature, p is the density of the fluid, p 
the viscosity, a is the thermal diffusivity, the coefficient 
of thermal expansion, g is the acceleration due to gravity, K is 
the permeability of the porous medium and A and B are 
respectively the width and height of the enclosure. 

The problem may now be restated in the dimensionless form. 
The dimensionless variables chosen are , 


A ’ y = A 


Tj. - T 

ft _ t s 

■ T - T 

b s 


u = 


UA 
a ’ 


VA PA^ ^ ^ 

= ar • p = 2 • ® = T^- r 

pa b s 


Pr = ^ 


pa 


. Da = ^ 


, As = 


B 

A ’ 


* _ - T^>KA 

Ra — 

pa 

* 


(3.6) 


where. Da and Ra are Darcy number and Ray 1 e i gh-Darcy number 
respectively , Pr is the Prandti number and As is the aspect 
ratio. In terms of these dimensionless variables, the equations 
<3.1> - (5.4) may be rewritten as. 


- 0 

dx dy 


(3.7) 


u 


du , ©u 

= - Jp + 

Pr( 

a^u 


ax 

ay 

ax 

ax^ 

ay^ 

dv , ©V 

+ V 

= - + 

Pr( 

a% 


ax 

ay 

ay 

ax^ 

ay^ 


) - 


Pr u 
Da 


(3.8) 


> ^ 

Da 


(Ra*© 


V) , (3.9) 
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aa ^ ae 

dx 


2 2 

a e as 

2^2 


<?. 10 > 


the boundary 

cond i t i ons . 


X = 0 : 

u 

= 0. 

II 

CD 

o 

II 

> 

0 . 

It 

X 

u 

= 0. 

V = 0 , e = 

0 . 

y = 0 : 

u 

= 0, 

V = 0 , e = 

1 . 

y = As : 

u 

= 0. 

11 

CD 

o 

II 

> 



( 5 . 11 ) 


Since the pressure field is unknown, it is eliminated from the 
equations <5.8> and <5. 9). Stream function ^ and vorticity co 
are defined as. 


aw 

u = ^ 

ay 


ay/ 


and 


a> = 


au 

ay 


aw 

ax 


(?. 12 > 




The equation of continuity (3.7) is automatically satisfied by 
(3.12) and the vorticity can be written in terms of y/ as. 


a w , a w 

o> = 


(5.14) 


ax ay 

The momentum equations (3.8) and (3.9) are transformed into the 
vorticity equat i on , 


2 2 

aw ao> ayt aco ,aco,aw. Pr, 

ay ax ax ay ^2 ^2 Da ax 

ax ay 

The energy equation (3.10) now takes the form, 
dy dx - dx ay = ^yZ ’ 


(5.15) 


<5.16) 



59 


The boundary conditions to be satisfied by y/ and & are, 
X = 0 . V = 0. . 0, e = 0 . 

X = 1 : V' = 0 . ^ = 0, © = 0 , 


- <?.17) 

y=0:w=0. J^=0, 0=1 . 

y = As ;y/=0. ^=0,0=0^ . 

The advantage of stream function vorticity formulation is 
that instead of working with the continuity equation and two 
momentum equations only two second order equations are to be 
considered and the number of unknown variables is also reduced 
by one with the elimination of the pressure term. The principal 
difficulty in this formulation is that, in general, vorticity is 
unknown apriori along solid boundaries and hence is not useful 
in situations with moving boundaries, for example in 
solidification or melting problems for which velocity boundary 
conditions must be imposed on the moving boundary. However, the 
problem under consideration may easily be treated with this 
formulation and it appears to require less computational time 
than the velocity pressure formulation described in Chapter IV. 


3.3 NUMERICAL PROCEDURE 


The governing equations 
differential equations for 
solution is natural. Finite 
using upwind differences f 
differences for other 


are coupled non-linear partial 
which an iterative procedure of 
difference equations are derived 
or convective terms and central 

upwind difference 


terms . 


The 
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approx i mat i on of the convective terms depends on the flow 
directions in which a derivative is replaced by a backward 
difference if the velocity in the concerned direction is 
positive and by a forward difference if that is negative. The 
resulting finite difference equation is only first order 
accurate but leads to a diagonally dominant matrix and to a 
numerically stable scheme independent of the step sizes. The 
governing equations <>.14> - (?.16> are approximated in the 
following difference forms. 


j - + ’"i + l.j ^ + V-i.i + l . 


Oi . 


<Ax>' 


(Ay)^ 


1 , J 

<5. 18) 


lAy/.i 

1 o> . . - . . , 

' , 1.3 l±l . 3 

[Av'. 

^ 2Ax 

2Ay 

^ Ax 

Pr( 

to . , . - 2co . . + 

1-1.3 1,3 

“L+l.j + 


(Ax) 



w . - <J> . ... 

< '•>. > 

Ay 


w. .+a>. ... 

1.3-1 1,3 1,3+1 


<Ay) 


^ < o,. . + Ra* + ) . 

Da 1,3 2Ax 


(5.19) 


2Ay Ax 


®i+l 3 

' J ) + _ * < 

2Ax 


e 


. . - e. ... 

1.3 I . 3 ± 1 . 

- Ay 


^i-1 . 3 3 ^ ^i+1 . j 

(Ax)^ 


e. . - - 2©. . + a. ... 

1.3-1 1.3 1.3+1 


(Ay)' 


(5.20) 


where, Av'. = .-v'-, , hyj . - rfj . ■, and Ax 

^1 ’^l+l.j "^1-1,3 ^1.3+1 ^1,3-1 

and Ay are the mesh sizes in the x and y directions 
respectively. The subscripts i and 3 correspond to x and y 


d i rect i ons . 
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The difference equations (5.T8) - (5.20) are solved by 
line-by-line method presented by Patankar C1980>. This method 
has the advantage of being simple and easy to apply, involving 
the solution of tridiagonal sets of equations along the lines 
parallel to x and y-axis alternately called horizontal and 
vertical sweeps. Starting from below, all the horizontal lines 
except those with prescribed boundary conditions are swept one 
by one during the horizontal sweep. For a fixed j, the unknown 
variables corresponding to the node points (i, j+1) and <i,j-l> 
are assumed to be known from an initial guess or previous 
iteration or the most recently available values. A tridiagonal 
system is set up for the unknowns corresponding to the nodes 
<i-l,j), <i,j> and <i+l,j> which is easily solved by tridiagonal 
algorithm. After the horizontal sweeps, vertical sweeps are 
performed starting from left to right in which the nodal 
variables at <i+l.j) and <i-l,j> are considered as known and a 
tridiagonal system is solved for the nodal values at (i,j-l), 
<i.j> and <i,j+l> for a fixed i. The process is repeated for 
several iterations till convergence at ail grid points is 
obtained. 


Assuming r^ 
expressed in a tr 
sweeps as follows 


2Ax 


and 


Ay 

i diagonal 


2 

form 


2Ay 

”Ax 

for 


equation <5.19> may be 
horizontal and vertical 
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Horizontal Sweep 

(- jAv^ J -Zr^Pr) a>._j . + <lAy/J + jAy/.j + 4rjPr + Ar^Pr 

Pr 

+ 2 Ax Ay -p—y o> . - 2r,Pr «■ ,, . 

Da 1 , j 2 1 +1 . j 

= 2 r, Pr ci> . , + ( lAy/. I + 2 r, Pr ) o>. . 

1 I . ) -1 I ^ 1 I 1 1 , J +1 

D o P F' 

- Ay ^ <© - e . > , f or Av' > 0 and Ay/ . > 0 . 

Da 1 + 1.] i-l,] ’^i 

(- 1 Ay/^ I -Zr^Pr) t^j _2 j ■*■ I + Ar^Pr 

+ 2 Ax Ay CO. . - 2r„Pr co. 

Da 1 , ) 2 1 +1 . J 

= (iAy/.! +2r,Pr)co. , +2r,Prco. 

1 i.)-i 1 1,3+1 

Do P 'P 

- Ay - .1 " ©• 1 •> . "for Ay*. > 0 and Ay'. < 0, 

Da 1+1 , J 1-1,3 3 1 

-Zr^Pr j + I ^r^Pr + Ar^Pr 

+ 2 S' “i.i * ‘ “i+i,i 

= ZrjPr ». j.j + < lAv-jl + ZrjPr) ..j 

Do P p 

- Ay _ , .) , for Ay/. < 0 and Ay/. > 0, 

Da 1+1 .3 1-1,3 J 1 

and 

-Zr-Pr . + <jAy/.| + |Ay/.| + Ar.Pr + Ar.Pr 

illl,J 1 3 ■! ^ 

Pr 

+ 2 Di' "i.i ^ “i+l.j 

= (jApJ + 2rjPr> + ZPjPr o>, 

Ra* Pr 

" Da ^®i+l,3 ■ ®i-1.3^ ' "^^3 ° ^ °- 


<5.21) 
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Vertical Sweep 

-ZrjPr ,_j + (jAv'j] + + 4rjPr + 4r2Pr 

+ 2 Ax iy g.) j + (-lAvJ - 2rjPr) 

= <|Av.j| + 2r2Pr) + 2r2Pr j 

Rp* Pr 

- Ay - © , > . for Ay/ >0 and Ayj . > 0. 

Da 3+1, j 1-1, j ’^i 

< - I A\y . j -2r j Pr ) w. ._j + < jAy'. | + [Av/.j + 4rjPr + 4r2Pr 

+ 2 Ax Ay £ 0 . - 2r,Pr 6> . ,, 

Da 1 , j 1 1,5+1 

= (jAv'.l + 2r-Pr) <;).,.+ 2r,Pr « 

'j* 2 1-1.} 2 1+1,] 

•it 

- Ay - e. , .) . for Av^. > 0 and Ay^ . < 0, 

Da 1 +1 , J 1 -1 , ] ^ ] 1 

- 2rjPr oj. + <jAv/.| + jAy' . | + 4rjPr + 4r2Pr 

* 2 Ax Ay g) <yj J + <-|A»-j| - 2rjPr) 

- Ay , . - 9. . .) , for Ay, < 0 and Ay, > 0, 

Da 1 +1 , J 1 -1 , J J I 

and 

<- j Ay . j -2rjPr> oj. + ( |Ay . j +• jAy^ j + 4rjPr + 4r2Pr 

+ 2 Ax Ay g) », J - 2rjPr 

= 'lAVjl + Zr^Pr) j + Zr^Pr ^ 


< 5 . 22 ) 
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Coefficients on the left hand side and ail the expressions on 

the right hand side in the equations <5.21) and <5.22) are 

evaluated from the previous iteration values of ^i-fl j 

and c*> , . The boundary conditions for o> are obtained using the 

1 "1 , } 

no-si ip condition at the right boundaries and the expression 


<5. 14) . 

At X = 0 : 


<jd 


Z 

\ 1 

dy. 




By 


= 0 and v' = 0 <no-slip). 


Applying central differences at a point <l,j), one gets , 
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2, j 




1 . j 


(Ax) 


■ ^O.i 
2Ax‘ 


0 and v'l • = 0 

1,1 


Thus 03 


1 . j 


(Ax)^ 


, f or al 1 j . 


Simi lariy , 


2V', 


at X = 1 ; fj = - , for all j, 

(Ax) 

2 

at y = 0 : o). - = ’? * ^ 

(Ay)"^ 


ZW 


i ,M 


at y = As: o>. , = ‘'’V • ail i 

(Ay)^ 


(3.25) 


where, the number of grid points in x and y directions are N+1 
and M+1 respectively. 

In a similar fashion, the equation (5.20) is rewritten for 


horizontal and vertical sweeps as follows. 
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Horizontal Swoop 


(-jAv'jl-Zr^) + jAv-.l + 4rj + 4r2> 0 . . 


®i+l. j 


2r 


for Ay/. > 0 and Ay'. > 0 , 


<-|Av/. j-Zr^) + lAv-j + 4rj + 4r2> 


e 


1 . ) 


^^2 ®i+l . j 


(jAy^J + 2 rj) 6 .^ ._j + Zr^ 6 . ^ 


for Aw. > 0 and Ay'. < 0, 
J 1 


+ <-jAy'.| - Zr^) + 

+ (lAy-J f 2rj> . 

for Ay/ ^ < 0 and Ay#. > 0. 


= 2r, e. . , 

1 1 , J -i 


and 


■ ^^z + i^^jl + ®i.i 

+ <-jAy'.j - Zr^) e.+i J 

= <|Ay'J + 2rj> + 2rj j+j . 

for Ay/j < 0 and Ay'. < 0. 

<5.24) 
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Vertical Sweep 


+ 4rj + 9 . . 

* <-|iV-i| - 2rj> 9. 

= < 1 Av 2 j + 2r > e , + 2r. 6 


i-1 . j 


2 i+1 , j • 


for Aw. > 0 and Ay/. > 0. 
J 1 


(-lAy/J-Zr^) ®j j_| + + 1-^V'j | + + At^} 


e . 


1 . } 


- 2r, e. . 


+ 2 r„ e. 


and 


1 ''i.j+l * 2=^2^ ®i-l.j ■ "‘2 ''i+i.j 

for Av'^ > 0, and Ay'. < 0, 

®i,i-l * ■* * “‘'z’ ®i,j 

* < - |Av-j| - 2rj) e, ,,j 

= 2r, + (|Av,,| t 2r,> 


for Aw. < 0 and Ay'. > 0, 

r J T'j 


<-|Av.J-2rj) e, j.j + + I^V'jl + ■‘■'j + Sj j 


for Ay'^ < 0 and Ay'. < 0 . 

<5.25) 

Again, assuming r = ^ , the equation <5.18) is subjected to the 
line by line procedure as follows: 



Horizontal Sweep 




_ 2 2 2 
= - r ¥'i i , - r <Ax> w . 

1 . J-1 . ) + l 1 . j 


<>.26> 


Vertical Sweep 


2 2 2 
r V'i : , - 2<l+r > v . + r v' .l, 
1 • ) i 1 . ) 1 . j + 1 


- '“i-l.l - '“i + l.j + “i.i 


<5. 27) 


The rate of heat transfer along each wail is determined 
from the temperature field. The local Nusseit number along the 
side wall is , 


Nu = 


ae 

ax 


<3. 28) 


x=0 


Along the bottom wall it is , 


Nu 


ae 

ay 



C3.29) 


and along the top wall, the local Nusseit number is defined as , 


Nu = . <5.}0> 
-y^As 

The derivatives on the right hand side of equations 
< 3 . 28 )-( 3 . 30 ) are computed by using a three-point finite 
difference formula. Once the local Nusseit number along the 
surfaces are determined, the average Nusseit number along each 
surface is evaluated by Simpson’s rule of integration. 



DO 


The iterative procedure adopted is summarized below, 

(a> The first iteration is started by guessing values for the 
fluid temperature 6, stream function y* and vorticity « at 
all the grid points. 

(b) Equations {3.24> and <3.25> are solved for the fluid 
temperature 6. 

(c> The vorticity w is obtained from equations (5.21) and 

<>.22>, using the & values from the previous step. 

<d> From the calculated values of in step <c>, the stream 

function yj is obtained from equations (5.26) and (5.27). 

(e> The local Nusselt number and the average Nusseit number are 
estimated from the fluid temperature field. 

(f> The new guess for &, y^ and w is evaluated from their 

current values using suitable relaxation parameter and 
steps (b-e> are repeated until convergence is reached. The 
convergence is presumed to have been achieved when 


: k+1 k ^ 

^i i ■ '^i j ■ -4 

i k+1 ; 

> ■ 

, I . J ; 

where, the superscripts k+1 and k indicate the iteration on 
which that value was obtained. 


3.4 RESULTS AND DISCUSSION 

The numerical solutions have been made on a HP 9000 
computer. To test the validity of the numerical procedure, the 
average Nusselt numbers along the side wall obtained by using 
the present method are compared with the results of Chen et 
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a 1.(1 987) in Table 1 . for Pr = 7.0 .As = 1.0 and lo"^ < Ra < 10^ 
in the absence of porous matrix , where Ra is the Rayleigh 
number. The results tabulated in Table 1 are obtained for 
dimensionless side wall temperature = 1.0, top wall 
temperature = 0.7055 and bottom wall temperature = 0.0. The 
agreement is satisfactory and substantiates the correctness of 
the computational work. 

At the second stage the boundary conditions for the present 
problem have been imposed by setting the side wall 
temperature = 0.0, bottom wall temperature = 1.0 and three 
different values of the top wall temperature = 0.0, 0.5, 1.0. 
Computations for different grid sizes were carried out and 
61 X 61. 81 X 81 and 101 x 101 grid fields showed insignificant 
change in the results. It was observed that a grid system of 
81 X 81 nodes is a good choice. The physical domain is limited 
by 0 < X < 1 and 0 < y < As . 

Again for the purpose of comparison, the same problem is 
solved by finite element method described in Chapter IV and the 
local Nusselt numbers along all the surfaces obtained by the two 
methods are presented in Tables 2, 5 and 4. The small 
difference between FDM and FEM is probably due to fine mesh size 
used in FDM calculations than that used in FEM. 

Flow patterns and isotherms for various parameters are 
shown in Figures 2 to 9. There is a symmetry in the vortex 
pattern about the midpiane of the enclosure, due to symmetric 
nature of the boundary conditions prescribed in the problem. 
The flow consists of a single cell filling the entire half of 
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the enclosure rotating in the clockwise direction but 
anticlockwise in the other half. The velocities are higher near 
the wall and lower near the line of symmetry. It is observed 
that the flow is almost stagnant near the corner. The 
circulation strength depends on Ra , Da, As and Art increase 
in Ray 1 e i gh-Darcy number increases the strength of the 
circulation vortex. It is seen that this change in the vortex 
strength is nearly proportional to the Ray 1 e i gh-Darcy number. 
As the Darcy number decreases, the strength of the vortex 
decreases due to increasing resistance of the porous matrix. in 
order to study the effect of aspect ratio, three different sizes 
of the enclosures are chosen. The patterns are very similar for 

As = 0.5, 1 and 2 except that the values of stream function and 

% 

temperature are larger for large aspect ratio. Streamlines and 
isotherms are also plotted for two different values of top wall 
temperature = 0.5 and 1.0. The change in the top wall 
temperature also alter the flow structure. The circulation 
pattern is smoother in the case of 6 ^ = 1.0 when compared to the 
case = 0.5. 

Furthermore, for small values of Ra, a pure conduction 
solution with parallel isotherms were obtained. As the Rayleigh- 
Darcy number increases, the effects of convection are noticed. 
For moderate values of Ra*. the isotherms are evenly spaced 
curves indicating that the heat transfer is dominated by 
conduction and the convection effects are less. The isotherm 
patterns are more distorted in the case of large Ray 1 e i gh— Darcy 
number when the convection is dominant. 



Figures 10 and 11 show the local Nusselt number along the 
bottom surface as a function of Ray 1 e i gh-Darcy number for 
different values of Da and ©^. The Nusselt numbers along each 
surface are found to increase with Ray 1 e i gh-Darcy number but 
decreases considerably as Da increases. Nusselt number along 
the bottom surface is the largest and it is the least along the 
top surface. Increasing increases the Nusselt number along 
the top surface. Increasing the aspect ratio, decreases the 
Nusselt number along the top and bottom surfaces. Along the 
side wall, increase in the aspect ratio brings a large increase 
in the Nusselt number for ©^ = 0. For larger magnitudes of ©^, 
changing the aspect ratio had little effect on the Nusselt 


number . 
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Table 1: Average Nxisselt nvimh^r ctloTtg th& sid& wclH in ai>s^nc^ 
of poraxis matrix for Pr = 7.0. As = 1.0. sidbe wall 
i&mp^ratxir& = 1.0, top wall temporalxiro = 0.705 5 and, bottom wall 
tomp&ratxjtro = 0.0. 


Ra 

Compar i son 


Left wall 

R i ght wal 1 

10^ 

Chen, Ho and 

Humphrey 

5.5409 

3.5408 


Present 


3.6005 

5.6005 

10^ 

Chen, Ho and 

Humphrey 

4.5244 

4.5245 


Present 


4.3506 

4.5507 

10^ 

Chen, Ho and 

Humphrey 

4.8472 

4.8472 


Present 


5.1757 

5.1757 


Table 2: Local Husselt nwnb&r along the side wall for 





Pr = 

1.0, As 

= 1.0, 

©t = 0.5 . 




Ra = 

10, Da = 

lO"^ 

s 

. Ra = 

10, Da 

= 10-^ j Ra* = 

100, Da 

= lO" 

y 

0.4 

0.5 

0.6 

1 0.4 

0.5 

0.6 ^ 0.4 

0.5 

0.6 

FDM 

1 . 539 

1.315 

1.258 

^ 1.514 

1.295 

1 

1.257 ; 1.895 

2.070 

2.279 

FEM 

1.551 

1.321 

1.260 

1.515 

1.305 

1.246 5 1.884 

2.092 

2.502 
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Table 5: Local Ni/ssoll nioribor along the top wall for 

Pr = 1 . 0 . As =1.0,©^= 0.5. 


Ra*= 10. Da = 10~^ | Ra*= 10. Da = lO'"^ ' Ra* = 100. Da = 10 



- 


1 



4 

- 


X 

0.5 

0.4 

0.5 ! 

0.5 

0.4 

0.5 0.5 

0.4 

0.5 

FDM 

0.919 

0.700 

0.640 ; 

0.925 

0.704 

0.640 ; 0.069 

0.058 

0.144 

FEM 

0.910 

0.710 

0.642 s 

0.916 

0.715 

0.644 i 0.064 

0.054 

0.158 


Table 4: 

Local 

Mussel t 

number 

along 

the bottom, wall for 





Pr 

= 

1.0, As 

= 1.0. 

e = 

0.5. 






-3 

[ 



-4 

i * 

100, Da 



Ra = 

10. Da = 

10 ^ 

V 

Ra = 10. 

Da = 

10 

; Ra = 

= lO" 

X 

0.5 

0.4 

0.5 

} 

0.5 0 

.4 

0.5 

i 

1 0.5 

0.4 

0.5 

FDM 

2.278 

1 . 841 

1.665 

i 

2.516 1 

.850 

1.714 

$ 5.500 

1.994 

1.515 

FEM 

2.290 

1.854 

1 .680 

> 

1 

2.516 1 

.861 

1.720 

i 5.610 

2.110 

1.551 
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Figure 1 The Physical system 








Figure 9 Streamlines and isotherms 
for e,=1,Pr = 1,As=1/2,R*a=10^ 

Da=10^ 




Figure 10 Local Nusselt number along the 
bottom wall for 6^ =1jPr =L As=1 




Figure 11 Local Nusselt number along the bottom 
wall for Pr =1 1 As=l > Da =10 



CHAPTER IV 


NATURAL CONVECTION IN A VERTICAL CYLINDRICAL ANNULUS 
FILLED WITH A POROUS MATERIAL 

4. 1 INTRODUCTION 

In spite of its great practical importance in many branches 
of engineering, the subject of natural convection in porous 
media with cylindrical geometries is not nearly as well studied 
as the rectangular counterpart. Very recently, free convective 
heat transfer in cylindrical annuli filled with fluid saturated 
porous media has been studied experimentally and theoretically. 
Reda <1983). Prasad and Kulacki (1985) and Prasad et 
al . <1985,1986) reported some experimental results on steady 
state thermal field throughout an annular region bounded by a 
vertical inner cylinder and a concentrically placed isothermal ly 
cooled outer cylinder for different thermal conditions on the 
heated inner wall. In a theoretical study, Bau and 
Torrance <1981) examined the onset of natural convection in a 
fluid saturated porous medium contained between vertical coaxial 
cylinders of isothermal horizontal boundaries with heating from 
below and cooling from above. Hickox and Gartiing (1982) 
considered natural convection in a vertical porous cylindrical 
annulus whose inner wall is heated at a constant temperature, 
the outer wall is isothermal ly cooled and the top and bottom are 
insulated. They obtained finite element numerical solutions and 
also used an approximate analytical method to obtain heat 
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transfer results for low Rayleigh numbers and high aspect 
ratios. Later Havstad and Burns <1982> analyzed the problem 
using three different methods, the finite difference numerical 
method, an approximate asymptotic scheme and perturbation 
method. Numerical results for heat transfer were presented for 
moderate cylinder spacings and high temperature difference. 
Perturbation results valid for all cylinder spacings at low 
temperature difference were obtained. Asymptotic solutions for 
very tall cylinders and all temperature difference were 
supplied. They also presented the effect of wall Biot number 
when the heat is rejected to the ambient through a conducting 
wall. The work of Philip (1982) is also applicable to natural 
convection in porous media at small Rayleigh numbers. Prasad 
and Kulacki (1984) reported numerical results for steady free 
convection in a vertical annulus filled with a saturated porous 
medium whose vertical walls are at constant temperatures and top 
and bottom are insulated. These authors considered the problem 
for a wide range of Rayleigh Darcy number , radius ratio 
parameter and aspect ratio ^ 1 . Curvature effects on 
temperature and velocity fields were shown to be significant. 
In another paper Prasad and Kulacki (1985) determined heat 
transfer rates in short cylindrical porous annulus both 
numerically and experimentally for the same thermal conditions 
on the boundaries. Recently Prasad (1986) obtained numerical 
results for the case when the inner wall is heated by applying a 
constant heat flux which results in a higher rate of heat 


transfer compared to the isothermal heating. 


In all these 
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theoretical studies Darcy’s law is assumed to hold and boundary 
and inertia effects have been totally neglected. 

Consequently, the present analysis deals with free 
convective flow in a vertical cylindrical annulus filled with a 
fluid saturated porous medium without neglecting these effects. 
The vertical walls are at constant temperature with the inner 
being hot and the outer being cold. The top and bottom of the 
annulus are insulated. The boundaries are impermeable and the 
results are based on the generalized momentum equation. The 
problem is analyzed via finite element method based on the 
Galerkin form following a procedure outlined by Taylor and 
Hughes <1981). 

4.2 THE BASIC EQUATIONS 


Assuming local thermal equilibrium between the solid and 
fluid, and also, a Boussinesq fluid, the axisymmetric 
two-dimensional form of the governing equations expressing 
conservation of mass, momentum and energy for the present 
problem are. 
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The geometry of the cylindrical porous annulus and the 
coordinate system (R,Y> are shown in Figures 1 and 2. U and V 
represents the velocity components in R and Y directions 
respectively, p is the density of the fluid, p the viscosity and 
ft is the thermal expansion coefficient of the fluid. K is the 
permeability of the isotropic porous medium and a is the thermal 
diffusivity. T, P and g denote the temperature, pressure and 
the gravitational acceleration respectively. Denoting the inner 
and outer cylinder radii by R. and R^ , the cylinder height by 
H, gap width of the porous annulus by D and temperature at the 
inner and outer walls by T. and ; the appropriate boundary 
conditions for the problem as indicated in Figure 2 are. 


R = R. : 

1 

U = 0. 

V = 0. 

T = 

T. 

1 

R = R : 
o 

U = 0. 

V = 0. 

T = 

T 

o 

Y = 0 : 

U = 0. 

V = 0. 

aj 

5y 

= 0 

Y = H ; 

U = 0. 

V = 0, 

ai 

dY 

= 0 


On introducing the dimensionless quantities. 


r = p , y = n ' ^ 


P 


2 

P« . 


, & = 


T - T 

o 

T.- T 
1 o 


R ^ pgftKDi T.- T„> ,, 

1 D ' D pa pa 


(4.6) 


the governing partial differential equations can be written in 
the following non-dimensional form. 
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The boundary conditions in the non-dimensional form are, 


r = r . : u = 0 , 

r = r.+l : u = 0 , 

y = 0 ; u = 0 , 

y = 0 ; u = 0 , 

In the dimensionless form, the radius ratio x = 1 + — , the 

r . 

1 

aspect ratio As , the Prandtl number Pr, the Rayleigh-Darcy 
number Ra and the Darcy number Da are the important parameters. 

4.3 DEVELOPMENT OF FINITE ELEMENTAL EQUATI(»IS 

First, the geometrically complex domain of the problem is 
represented as a collection of geometrically simple sub domains, 
called finite elements . Associated with each element are eight 
discrete node points, each of which is specified by a local 
co-ordinate as well as a global co-ordinate system. The eight 
noded isoparametric elements are chosen to accommodate parabolic 


V = 0 , © = 1 , 


V = 0 , © = 0 , 


V = 0 , ^ = 0 , 

V = 0 , ^ = 0 . 

dy 


<4. 11) 
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variation in the variables on the boundary and in the interior 
of the element. In such elements, the temperature and velocity 
components are specified at all the eight nodes whereas the 
pressure is specified only at four corner nodes. The shape 
function (also called as interpolation functions) for eight 
noded element in terms of local co-ordinates ? and n are given 
by the following formulae. 

For corner nodes ; 

Ni = j (1 + + r>.T 7 )<f + 77.77 - 1) . i = 1.5. 5. 7 . 

<4.12a> 

For midside nodes : 



1 - 1 

+ 77j77>,if = 

0 . 

i = 

2.6 

<4, 

.12b) 

N. = ^ ( 

1 +?.?>( 1 

- 77 ^), if 7?. = 

0 . 

i = 

4,8 



For pressure, linear 

interpolat i on 

is used 

and 

the 

shape 

functions corresponding 

to four corner 

nodes 

are 

* 



Ml = i < 

1 - >( 1 

- 77 > . 






- 5 ' 

1 + f )< 1 

- 7? ) . 
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> 4 

M. = 4- < 1 - ? >( 1 + 77 ) . 

4 4 ^ 

The advantages of using eight noded isoparametric elements and 
depicting the variation in pressure by shape functions of one 
order lower than those for defining the velocity distribution 
are discussed in detail by Taylor and Hughes (1981). 
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In terms of these shape functions. the velocity, 
temperature and pressure over an element are written in the 
following form, 

8 

u = E N u . 
i=l * * 

8 

V = E N. V. . 

itl ‘ * 

8 <4.14) 

e = z N 0 

i = l ' ‘ 

4 

P = E M p . 
i=l * * 

Spatial coordinates within the element are defined in terms of 
all eight nodes. Using these quadratic and linear approx irrat ion 
functions, and employing Galerkin weighted residual approach in 
which the weighting functions are taken to be the shape 
functions themselves, equat i ons (4 . 7 > - <4. 10) become, 

ne dN. N. aN 

r rr M.< « ’ u. + -1 u . + - ’ V. >r dr dy = 0 . <4. 15a) 

" JJ 1 ar j r j ay j 
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dN 
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The repeated suffixes indicate summation with the suffix running 
through the integral numbers 1,2, such that 


and 


I = 1. 2, 5, 4 , 

j = 1,2, .8 

k = 1,2, ,8 . 


The outer summation refers to each element in the domain. The 
symbols u^ , v. , 0. and etc. denote the nodal values of u, 
V, O and p respectively within a particular element. 

Invoking Green’s theorem the second order terms in the 
equations (4.15b> - <4.15d> are reduced and the finite elemental 
equations for mass, momentum and energy balance are written as 

ne dN . N . dN . 

E // ”i‘ af’ “j ’ “j ’ ay ''i = o , (A.Ua) 
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r 

where r denotes the element boundary. In the nonlinear 

inertial terms IT^ and represent the velocity components at 

the kth node which are taken to be known (from the previous 
Iteration) in order to make the equation quasilinear. The 
boundary integrals on the RHS of equations <4.16b> involve the 
derivatives of the velocity components and the temperature in 
the direction of the outward normal to the boundary. However. 
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for the present problem with boundary conditions given by 
(4.11). the RHS integral contributions are zero. 

4.4 MATRIX FORMULATICaf AND S(X.UT1CM PROCEIXIRE 

The finite element equations <4. 16a) - (4.16d) are 
assembled in the global matrix equations to solve the nodal 
values of pressure, velocity components and temperature. The 
assembled matrix equations take the form. 


A X 
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t 
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where, the chosen form for 
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and A is generally denoted 

as 

the 

stiffness matrix in 

the 

FEM 

terminology. Due to the presence 

of the inertial terms 

i n 

the 

momentum and energy equations 

, A is not symmetric. Because 

of 

the nonsymmetry of this 

matrix. 

the memory space 

and 

the 

inversion time are larger. 

Also 

the inertial terms 

i n 

the 


governing equations make the matrix equation (4.17) nonlinear, 
hence an iterative solution procedure is required to be imposed. 
Using u^ and v^ values from the previous iteration. the 
stiffness matrix A is updated and the solution vector X which 
contains all the nodal values of u, p, v and 6 is calculated. 
This is continued till the results do not change significantly 

from one iteration to the next. 

The solution procedure used here for solving the assembled 
matrix equations is known as the FRONTAL METHOD . This method 
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leads to a significant reduction in the con^uter storage 
requirement and the execution time. In the present problem the 
assembled elemental coefficient matrix is itself large with a 
size of 28 x 28. Also there is a limitation on the minimum 
number of elements to be used to discretize the domain, in order 
to avoid the i 1 1 ~condi t i on ing of the global assembled matrix 
during inversion. A fairly large number of elements are used. 

The frontal technique originally devised by Irons <1970). 
has proved to be very effective for solving positive definite 
matrices. This method has been modified by Taylor and 
Hughes<1981> for solving steady state Navier Stokes' equations. 
Originally this method was based on direct Gaussian elimination 
for solving symmetric matrices, where the leading diagonal is 
always used as pivot. For unsymmetric matrices encountered in 
wide range of problems in engineering fluid mechanics, the most 
suitable pivot is not always the leading diagonal, and frontal 
solutions exist for off-diagonal pivoting. These, however, tend 
to be more time consuming. The method used here uses only 
diagonal pivoting incorporating many features of the frontal 
method for solving symmetric matrices. 

In Frontal method the reduction is carried out on an 
elemental basis, to avoid the difficulty of large storing space. 
The complete global matrix is never assembled fully. As soon as 
the contributions from all the elements to a particular nodal 
point have been completed, the corresponding variables 
associated with the node are eliminated. A preassigned global 
matrix core area is filled from contributing elements, the 
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largest diagonal entry in this is then found and used as a pivot 
in the direct Gaussian elimination process. The corresponding 
reduced equations written onto disk and more elements and 
corresponding equations are taken into the core. When all the 
elements have been included and the variables reduced, the back 
substitution process is initiated (Taylor and Hughes 1981). 

The solution procedure involves the following steps: 

(i) Calculate the elemental coefficient rratrix CA]--, -- and 

7L oxZ o 

For the element under consideration. 

<ii) Assemble the matrices Into the global coefficient matrix 
and global right hand side vector. 

<iii) Introduce boundary conditions if any, for the particular 
e 1 ement . 

<iv> Reduce the global matrix using the Gaussian elimination 
procedure . 

<v) Repeat the steps <i) to <iv) for the next element. 

4.5 RESULTS AND DISCUSSION 

Numerical results using finite element method have been 
obtained for a wide range of dimensionless parameters. The 
solution domain has been discretized by using 100 rectangular 
elements, 10 in each direction. To test the applicability and 
accuracy of the present method, the results were obtained for 
the vertical cavity case and were compared with the solution 
obtained by finite difference method described in Chapter III. 
Both the methods yield comparable results confirming the 
validity of the present scheme. 
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The governing equations of the present problem are solved 
to first obtain the velocity and temperature fields. 
Subsequently the local Nusselt number &nd the average Nusselt 
number along the heated wall are evaluated using the expressions 
w ^ I 

y “ dr lr=0 ’ <4,19> 

and 

~ Xq Wr ^^lr=0 ' <4. 20) 


respectively. These are evaluated by the method of numerical 
differentiation and integration. 

The velocity and temperature fields for some typical values 
of Ra , Da, As and h are presented in Figures 5 to 7. It is 
observed that the magnitude of u velocity (in the radial 
direction) first increases with distance away from the heated 
wall and then it decreases to zero. This feature can be 
explained in terms of the existence of a velocity boundary layer 
adjacent to the hot wall due to the buoyancy effects and the 
entrainment of the enclosure fluid into the boundary layer. The 
entrainment velocity will obviously have a maximum value at the 
edge of the boundary layer. The U“velocity is negative in the 
lower half of the enclosure and positive in the upper half. 
This indicates that the flow approaches towards the hot surface 
at the bottom of the enclosure and the hot fluid moves away from 
this surface at the top. This sense of flow circulation is also 
confirmed by velocity vector plots (Figures 10-12). 
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In Figure 4, it is seen that the magnitude of the velocity 
at the same location increases with radius ratio *e. The reason 
for this is that the ratio of area between the outer and inner 
cylinders increases with «. Since at steady state same amount 
of heat flows from the hot wail to the cold, the heat flux on 
the hot wall is much larger < since its area is smaller > than 
that of the cold wall. The higher heat flux causes more 
buoyancy and it is therefore obvious that the velocity should 
increase with radius ratio. This fact is illustrated by the 
velocity vector plots ( Figures 10-12). 

The temperature distributions in Figures 5 to 7 clearly 
show the existence of a thermal boundary layer near the hot 
wall. Since the radius ratio is larger than one for all these 
plots, the heat flux near the hot wall is much larger than at 
the cold wall, thereby establishing a boundary layer type of 
temperature variation there. 


It is seen in Figure 5, that the thermal boundary layer 
becomes steeper as Ra is increased. This confirms the known 
trend that as the effect of buoyancy increases (higher Ra) the 
convective effects become stronger and the boundary layer 
becomes thinner. In the vertical direction, the boundary layer 
thickness increases with distance y as expected. The plots 
indicate that apart from steep gradients near the hot wail, 
moderately high gradients exist near the cold wall also. This 
is because the circulatory nature of the flow causes the hot 
fluid adjacent to the inner cylinder to come into contact with 
the cold wall at a subsequent time. It is also observed that 
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for radial positions which lie in the middle of the cavity the 
temperature variation in the radial direction is negligible. 

As the radius ratio is increased (Figure 6). the boundary 
layer becomes steeper adjacent to the hot wall and 
simultaneously, the gradients at the cold wall decreases. Both 
the trends are the consequences of the higher area ratio between 
cylindrical surfaces at larger «. with respect to the aspect 
ratio (Figure 7). however, it is observed that boundary layer 
becomes thicker, when As is increased. The temperature 
gradients on the cold wall also decrease with aspect ratio. 

Local Nusselt number variation along the hot wail is 
plotted in Figures 8 and 9. The Nusselt number decreases with 
vertical distance since the thermal boundary layer thickness 
increases with distance. 

For each radius ratio and Darcy number it is seen that 
Nusselt number increases with Rayleigh number. As explained 
earlier, this occurs due to the boundary layer thickness 
decreasing with Rayleigh number. It is interesting to note that 
the maximum Nusselt number does not occur at y = 0 (bottom 
corner) but some distance above it. This is believed to be a 


consequence of the existence of the stagnation region in the 
corner zone with a decrease in the value of ae. The Nusselt 

number decreases obviously due to smaller ratio of the cold wall 

M- ■ 

area and hot wall area. Final ly, ^with increase in Darcy number, 
the Nusselt number decreases. This is because the viscous and 

Ife-' i ■ ^ . 

thermal boundary layer thicknesse# decrease when Darcy number is 
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decreased (higher resistance offered by the solid matrix). A 
thinner boundary layer, in turn, leads to larger Nusselt number. 

The velocity vector plots for various parameters are shown 
in Figures 10 to 12. It is observed that the velocity values 
adjacent to the hot wall are always larger than those near the 
cold or adjacent to horizontal sides. This is the result of 
larger buoyancy force existing near the hot wall which has a 
smaller area near the cold wall. The magnitude of velocity 
decreases as one approaches the mid portion of the cavity. It 
is not surprising. since the natural convective flow is a 
phenomenon which is generated by buoyancy effects at the 
boundary walls and in the interior portion of the cavity the 
buoyancy effect is weak. The strength of the flow is observed 
to increase, as expected, with the radius ratio. Also, the 
magnitude of the velocity vector is larger, especially adjacent 
to the hot wall. Clearly this is the consequence of larger heat 
flux arising due to smaller heat transfer area, with the 
decrease in the value of Jt. If the aspect ratio of the cavity is 
increased, then also the strength of the natural convective flow 
increases since the contact length with the non-adiabatic 

surfaces and the fluid increases. 

Numerical results for Darcy law may be recovered by setting 
Da very small. The magnitude of the Darcy number in the present 
study is in the range < DaSlo’'. Table 1 presents results 
of a parametric study of average Nusselt number along the heated 
wall as a function of Ra*. Da. . and As. Results are compared 
with those obtained for a Darcy flow regime by Prasad and 
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Kulscki (1984) • An incrG3i66 in 4 ‘Ko d i * i_ 

.icrease m the Ray 1 e 1 gh-Darcy number 

increases the average Nusselt number and an increase in the 

Darcy number produces a reverse trend. As the Darcy number 
increases, present results deviate from the results presented by 
Prasad and Kulacki (1984). The difference between the present 
results and the results corresponding to the Darcy flow is due 
to the inclusion of viscous and inertia terms in the momentum 
equation. Vafai and Tien <1981) have pointed out that the 
Nusselt number decreases when the inertia and boundary effects 
become significant. Results presented in Table 1 supplements 
this observation. 


Table 1: Present Husselt ntwvbers compar&d with the results of 
Prasad and Kulacki 


« = 2.00 ae = 9,0 


Rs 

90 

79 

100 

100 

80 

120 

As 

1 

1 

2 

2 

2.9 

2.9 

Prasad and Kulacki 

2.77 

5.99 

4.09 

9.72 

4.83 

9.79 

Present results 







with Da = 10~^ 

2.71 

3.91 

3.94 

9.60 

4.79 

9.67 

Present results 
with Da = lO"^ 

2.69 

5.39 

3.83 

9.49 

4.68 

9.96 
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Figure 11 Velocily Vecior for Pr = 1*0, Ra=100, Da =10, As=1-0 

(1mm length of arrow =1 unit of velocity) 




Figure 12 Velocity Vector for Pr =1*0 , Ra =250, Da =10 » K=4*0 

(Imm length of arrow =1 unit of velocity) 



CHAPTER V 


A NUMERICAL STUDY OF OSCILLATORY CONVECTIS^ FLOW 
THROUGH A POROUS MEDIUM 

5.1 INTRODUCTION 

The problems of steady free or mixed convection adjacent to 
vertical fiat plate embedded in a porous medium have been 
investigated extensively. Little attention, however, has been 
given to transient convection in a porous medium. The first 
paper on transient boundary layer flow in a porous medium is due 
to Johnson and Cheng (1978) who have obtained similarity 
solutions for specific variations of wall temperature in both 
time and position. Cheng and Pop (1984) presented approximate 
solutions for the transient free convection boundary layer in a 
porous medium adjacent to a vertical semi -inf ini te plate with a 
step increase in wall temperature or heat flux. Raptis (1985) 
considered two-dimensional unsteady free convective flow through 
a porous medium bounded by an infinite vertical plate, when the 
temperature of the plate is oscillating with time and obtained 
solutions with the restriction that the amplitude of oscillation 
is very small. Singh et al . (1986) studied the same problem 
without any restriction on amplitude by developing two 
asymptotic expansions in powers of the frequency parameter using 
a regular expansion technique for small frequency parameter and 
the method of matched asymptotic expansion for large frequency 


parameter . 
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Recently Raptis and Perdikis (1985) analyzed the effect of 
free convection on How through a highly porous medium bounded 
by an infinite ver cal wall of constant temperature absorbing 
the fluid with constant velocity when the free stream vibrates 
about a constant mean. The authors have shown their interest in 
the problem where the amplitude of free stream oscillations are 
very small. Nevertheless we may come across physical situations 
where the amplitude of oscillations will not be necessarily 
sma 1 1 . 

With this motivation, the present work is devoted to the 
study of convective flow and heat transfer along an infinite 
vertical wall of constant temperature. The free stream 
temperature is also constant but differs from that of the wall. 
Due to differences of temperature, density gradients are setup 
giving rise to buoyancy forces. The fluid is being constantly 
absorbed by the vertical wall. The steady motion is disturbed 
due to oscillations in the free stream at t = 0. The velocity 
field subsequently differs from the steady state solution. The 
free stream velocity has been chosen of the form, 

u = U (l+£:sin&)t>, 
o 

with no restriction on e except that it must be less than 1. 
The problem is solved by finite difference method using Crank 
Nicolson scheme. The flow behaviour has been computed for 
various values of s- Knowing the numerical values of the 
velocity field, the skin friction at the wall is computed using 
Newton's interpolation formula. The results obtained by this 
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method for small c agree satisfactorily with those given by 
Raptis and Perdikis (1985). 


5.2 GOVERNING EQUATIONS 

Consider an infinite vertical flat plate, with X-axis along 
the plate and Y-axis normal to it. The vertical plate is at 
constant temperature and absorbs the fluid at a constant 
velocity while the free stream velocity vibrates about a 
constant mean. We assume that all the fluid properties are 
constant except the density variation with temperature in the 
body force term. the governing equations for the problem are 


(Raptis and Perdikis 1985). 

^ = 0 


(5.1) 


au . au *^*^00 ^ ^ . 

M' ^ ay ~ dt' 9 " ^00^ 




o ai ^ ^ a^T 

^ Wt' ^ av " ' 


(5.2) 

(5.5) 


The boundary conditions of the problem are. 

Y = 0:U=0.V=-Vq= constant. T = T^. 

Y -► 00 : U = (1 + sin «' t' ) . T T^ . (5.4) 

Here. (U.V) are the velocity components along the axes, g is the 
acceleration due to gravity, ft is the coefficient of thermal 
expansion, v the kinematic viscosity. K' the permeability of the 
porous medium. Uq a constant velocity. S is the heat capacity 
ratio of the fluid filled porous medium to that of fluid. « the 
thermal diffusivity. T is the fluid temperature. T^ and T^ 
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the surface temperature and temperature of the free stream and 
and w' are amplitude and frequency of vibration of the free 
stream velocity respectively. 

Introducing the dimensionless quantities. 


U, 


u 


0 

V(X>' 


U„ T - T 

® e = - 00 . 

• j' _ j > t 

0 a 


t' V, 


Y V. 


U 


. y = 


I'g/? <T» - T ) Vn 

6> = -^ . Pr = ~ , Gr = -. “ - , K = ° K' 

'^0 “o % 


and substituting for the solution of the continuity equation, 
V = ~ '^0 ~ constant , 

equations <5.Z) and (5.3) become, 

* 


du £u _ du’^ d u . 1_ y * 

at " dy " dt 'y2 K 


u> + Gr © , 


d© _ a© _ 1 d^e 


FF . 2 • 

dy 


at ay 

The boundary conditions reduce to. 


<5. 5) 


<5, 6) 


y = O:U=O,0=1, 

y -» co:u-»l+«sin(i)t, ©-fO 


(5.7) 


5.3 SOLUTION OF THE TEMPERATURE FIELD 

Since the suction velocity is constant, the temperature 
field remains unperturbed and hence the equation <5. 6) takes the 

form. 


„ a© ^ a^© _ n 

Pr 3- + , = 0 

ay 2 


<5.8> 


The solution of (3.8) together with the boundary conditions 


(3.7) is g i ven by . 

e = .-’’P >' 


(5.9) 
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g, 4 numerical SOLUTION FOR THE VELOCITY FIELD 

Substituting the expression for © from <5.9>. the equation 
(>.5) is reduced to. 

du ©u _ du* , 1 , * -Pr y 

at ay ~ dt : 2 - u> + Cr e . (5.10) 

dy 

Equation (5.10) is to be solved numerically, together with 
the boundary conditions, 
y = 0 : u = 0 , 

y->oo;u->l+tfsina>t. (5.11) 


The initial condition is obtained by solving the steady 
state equation corresponding to (5.10) which is given by, 

. . , . ' Gr ’’l’' 

t . 0 : “ • 1 + ; <p,, , <pr , pp - ' 


■Pr y 


Gr e 


(Pr + Rj) (Pr + Rg) ’ 


(5.12) 


where, R, 


-1-|1 + 47k 


and R, 


•1+il + 4/l< 
2 


‘1 2 2 
The finite difference representation chosen for equation 
(5.10) is. 


i , j 1 " ^i + 1 . j+1 ~ . j+1 . j • j \ 

At ■ ' ~ Z\ Ay Ay f 


u . ... - u , 

1 , j + 1 


£ to COS <0 jAt * j- * 

^ V 

^i-l.j i.) 1+1. J ^ 


1 " ^i-l.j+1 " ^^i.j + 1 ^i + 1.3+1 

2 


(Ay)‘ 


(Ay)' 


1 ^ 

+ \ 1 + ^^ sin to jAt 


“i.i.i * “i.) \ , 
2 


-Pr iAy 


Gr e 


(5.15) 
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and the initial and boundary conditions take the form. 


u ; 


i .0 


= 1 + ■< 


Gr 


(Pr + Rj> (Pr + R^) 


- 1 t e 

J 

-Pr iAy 


Rj iAy 


Gr e 

'<Pr + Rj> (Pr + R^) 


Urt . 0 , 

0 . J 


u . = 1 + £ sin » jAt 
n+1 . j 


(5. 14) 


The subscripts i and j appearing in (5.15) and (5.14) correspond 
to y and t respectively. 

Equation <5.1?> written for i = Hl>n constitutes n linear 
equations in n unknowns, u. where n refers to the number of 

rnesh points on y-axis. These equations may be written in the 
matrix form in which the coefficient matrix is tridiagonal and 
the method of Gaussian elimination is used to obtain a solution 

""i. j + 1- 

Knowing the numerical values of velocity field, the ski ^ 
friction at the wall is calculated. In dimensionless form it is 


given by , 


T = 


du 


y = 0 


This has been calculated by hu^erical diffsrehtiatich usihg 
Newton’s interpolation formula. 

5.5 RESULTS AND DISCUSSION 

, .ird.1 results obtained by Rapt is and 
The approximate analytical 

11 values of ^ show that u vary 
Pardikis (1985), for very small values 

7 5 Hence any value of y after 
insignificantly after y • • 

H to y computation time. 

may be chosen to correspon 
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8.0 has been chosen to correspond to y ^ in the entire 

computation. Ay and At are taken to be equal to 0.1 and 0.005 
respectively. The solution procedure is repeated to get 
velocity field at different time step until t = >.0. The 

computations are carried out for Pr = 0 . 71 ; Gr = 1 . 0 , 5 . 0 ■ K = 

1.0, 4.0, 8.0; o> = 5.0 and e = 0.1, 0.2, 0.5, 0.8 on a DEC 1090 
computer . 

To see the correctness of taking y = 8.0 corresponding to y 
00 , we have also run the program for other values of y such as 
8.5 and 9.0 corresponding to y oo and we get the same numerical 
values for u. To judge the convergence of the numerical method 
used, different values of At such as At = 0.002, 0.005 are 

chosen and no significant changes in the results are observed. 

For the purpose of discussing the results, the velocity 
profiles are shown in Figures 1 and 2. The velocity field 
increases rapidly, attains its maximum and then decreases and at 


y 2 k 8.0 there is no significant variation. 

In Figure 2, the velocity profiles are drawn for Pr = 0.71, 
Gr = 1.0, K = 4,0, to = 5.0, t = 1.5 and for four different 
values of c = 0.1, 0.2. 0.5. 0.8. The velocity profiles have 
the same shape for all values of e but the maximum velocity 
which corresponds to y 2* 2.0 increases significantly with the 


increase in e. 

Figure 5. shows the variation in velocity with y for 
different values of Br and K, when Pr - 0.71, ^ 0.8, w 5. 

and t = 1.5. The velocity is shown to increase with Gr and K 
and the difference of the velocity are greater as Gr increases. 
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Figure 4 illustrates the variation in skin friction at the 
wall with time. Curves are drawn for Pr = 0.71, Gr = 10 K = 

4.0, « = 5.0 and « = 0.2. 0.5 and 0.8. The skin friction at the 
wall is oscillatory in nature and the amplitude of oscillations 
increases with e. 

The results obtained by the present numerical technique for 
small values of e are compared with the analytical approximate 
solutions of Raptis and Perdikis (1985) in Table 1. All the 
numerical results are calculated for the case in which the free 
stream velocity is of the form. 

LJ = U (1 tcsinwt) . 

00 o 

The values of velocity distribution are found within 5 percent 
difference when c - 0.1 but the difference increases slightly 
when e = 0.2 as can be seen from Table 1. It is evident that as 
e increases the error involved in the analytical approximate 
method is more. Therefore, the present numerical method is 
preferable to the analytical method for solving this problem 
even for very small values of e. There is no solution available 
in the literature to compare present numerical results for 


larger values of c. 
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Table 1: Comparison of 

velocity dtciTthMlton ohtatn&d. 

by the 

approximate analytical 

method and the 

present numerical 

TUBthod 

for Pr = 0.71 . 

Gr = 1.0, K = 

1.0, t = j . 


- 

u 



s y 

Analytical 

% difference 

Numerical 

0.2 

0.4256 

0.4111 

2.95 

0.5 

0.8187 

0.8004 

2.24 

1.0 

1.1128 

1.0982 

1.51 

1.5 

1.1960 

1.1876 

0.70 

0.1 2.0 

1.2005 

1.1965 

0.55 

5.0 

1.1554 

1.1550 

0.05 

4.0 

1.1164 

1.1167 

0.05 

5.0 

1.0940 

1.0942 

0.02 

6.0 

1.0825 

1.0825 

0.00 

0.2 

0.4512 

0.4571 

5.15 

0.5 

0.8705 

0.8494 

2.42 

1 .0 

1.1807 

1.1651 

1.49 

1.5 

1.2677 

1.2570 

0.84 

0.2 2.0 

1.2725 

1.2670 

0.45 

5.0 

1.2266 

1.2258 

0.06 

4.0 

1.1875 

1.1874 

0.01 

5.0 

1.1647 

1.1650 

0.05 

6.0 

1.1551 

1.1551 

0.00 
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CHAPTER VI 


APPLICATION OF A THERMODYNAMIC METHOD TO BENARD 
CONVECTION IN HORIZONTAL POROUS LAYER 

6.1 INTRODUCTION 

Present chapter deals with the breakdown of stability of a 
layer of fluid subject to a vertical temperature gradient in a 
porous medium bounded by two horizontal rigid walls. The Benard 
convection is the name used to describe the motions arising from 
the thermal instability of a thin horizontal layer of fluid when 
a steady temperature gradient is maintained across it. When a 
fluid filled porous layer is heated from below, a somewhat 
similar physical problem arises and was first studied 
theoretically by Lapwood <1948). Reviews of the contemporary 
investigations of the onset of convection in porous media have 
been presented by Cheng (1978), Nield (1985) and Bejan (1987). 
Recently. Lebon and Cloot (1986) employed extended irreversible 
thermodynamics to study natural convection in a thin porous 
layer heated from below. With few exceptions, most of the 
studies ar« besed on either Darcy's law or Brinkman model. 

The main aim of the present study is to apply a 
thermodynamic method to obtain neutral stability curve when a 
steady temperature contrast is maintained across a thin 
horizontal saturated porous medium bounded by two rigid wails. 
The generalized momentum equation is used to analyze the 
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problem, which has no exact solutinn . . . 

solution and a variational technique 

based on the governing principle of dissipative processes (GPDP> 

proposed by Gyarmati (1969. 1970) is employed. The principle is 

quite applicable to various dissipative systems and is widely 

used to study the flow and heat transfer by Vincze (1975). 

Stark <1974). Singh <1976. 1980) and Singh and Raj (1985). The 

most general form of GPDP is represented by 

6j<cr-v'-<^>>dV = 0, (6.1) 

V 

for any instant of time under the constraints that the balance 
equations , 

p. + VJ. = .T,, <i = 1.2 f) , <6.2> 

are satisfied. Here cr denotes the entropy production per unit 
volume and unit time. yt and are the local dissipation 
potentials and the integration’ is considered over the total 
volume V of the continuum, p. is the partial time derivative of 
the density p.^, 3. is the corresponding current density and o-. 
is the production per unit volume and unit time of the various 
attributes of the system. The entropy production, f, in the 

case of irreversible processes taking place in continuum can 
always be written in the following bilinear form , 
f 

a=rj-X.>0. (6.5) 

i = l ‘ ‘ 

where. J. and X. are the thermodynamic currents and forces 
respectively. Due to the second law of thermodynamics cr is a 
positive definite quantity. According to the Onsager’s linear 
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theory of irreversible processes, the currents J. end the 
dissipative forces X. are given by the following constitutive 
laws , 

f f 

* k5l ^ik -^k ’ ^-2 •f> • 

< 6 . 4 ) 

where, the phenomenological coefficients L., and R. represent 

Ik ik 

conductivities and resistances, respectively and are constants. 
They satisfy the famous Onsager's reciprocal relations , 


and 


*-ik " Si • '^ik ' Si • 

their matrices are mutually 

f f 

^im Sk “ ^ *^im *"mk 

m-1 m=l 


< i . k = 1.2 f > , 

reciprocal . i .e. , 

<5ik . <i,k = l,2 


(6.5> 


f > . 


< 6 . 6 > 

Here, is the Kronecker symbol. The local dissipation 
potentials are the homogeneous quadratic functions of forces and 
currents respectively. These are defined as (Gyarmati 1970), 


V = T i j^i ^K^i-x. ^ ° • 

' I . i, "iK 

1 f IS *“ i 

Since in the case of transport processes the forces can always 
be generated as the gradients of certain scalar parameters f'j » 
the thermodynamic forces, may be expressed as . 


X. = vr. . 

1 1 


< 6 . 9 ) 
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Using (6.5). (6.7). (6.8) and (6.9) the principle 
the form , 


( 6 . 1 ) 


takes 



f 

E J • tt. 

= 1 » I 


f 

E 

i .k=l 


ik 


Vr.-TT. 


1 

2 


f 

E 

i ,k=l 


'ik 


dV 


0 


( 6 . 10 ) 


6. a THE DUAL FIELD METHOD 

The two sets of independent variables CJ^. and 

are connected with each other by the 
relations (6.4). In the dual field method we assume one set of 
these variables and obtain the other set with the help of the 
constitutive relations. In the irreversible transport 

phenomena, the variables F. are fundamental ones, since their 
gradients, driving forces of dissipative transport 

processes. We, therefore, approx irrate the set CVT^ , TT^ ^f^ 

» # » 

by another set CVTj , '^ 2 * • • • -. VT^]. Consequently the 

corresponding current densities are obtained by the following 
constitutive equations , 


J . = 

1 


f 

E 

k=l 


L. 


ik 


CFj, r^. 


K 


(i = 1,2 f) 


( 6 . 11 ) 

It is remarkable that the duality property of the governing 

principle is preserved and the two sets of fundamental variables 

r. and r* coincide in the case of exact solution. The principle 
1 1 

(6.10) with the help of (6.11) takes the form , 
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< 5 ; 


V 


i . k=l 


2 Lik<'^i-vr*>-<7r.-7r*> dv = o . 


( 6 . 12 > 


which together with the balance equations 

P- <r,> + vr*) = <. <r,>. <i = i.z f, , 

( 6 . 1 ?) 

serves the basis for the dual field method. We can. again, 
approximate the fundamental field in the following form . 


,<n> 


n 

= a.<t) g. <x> . 


.n 


<6.14) 


where, Cg. <x>]j_j are a set of linearly independent functions 
which satisfy the boundary conditions imposed on F. The 
coefficients, a.<t> are the variational parameters to be 
determined with the help of the principle. The balance equation 
takes the form , 


+ v.<L(r‘">> vr*'"') = <r'">) . 

wt 


(6.15) 

where, the current density J is replaced by the constitutive law 


J = L<r) VT , 


(6.16) 


which, follows from the general constitutive laws (6.11). The 
balance equation (6.15) serves the purpose to determine the 
field r* in the form and the solution of (6.15) with the 
appropriate boundary conditions may be written as , 


The 

the 


r* B r*<"’ = r*'"’ (i, t 

volume integral (6.12) i 
real physical processes. 


s maximum at any 
that is; for the 


. . . . a ) . (6.17) 

n 

instant of time for 
exact value of the 
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parameters P. and the curronf 

j ..u tne current densities J.. it is fundamentally 

important that the maximum is zero for any time (Gyarmat i 1969). 

In the application of approximate procedure, the volume integral 

generally becomes a function of time and therefore the volume 

integral may be integrated over the time interval 0 < t < oo 

during which the process is considered. Thus the principle 

(6.12) becomes , 

00 


<5 J J " J L - 7r*)» (TP - 7P*>' dV dt = 0 , 


(6.18) 


0 V 


In (6.18) we have confined our treatment for the case of one, P. 
parameter. The total variation of the principle (6.18) with the 
values of P and P* from (6.14) and (6.17) becomes , 


J 


" a?*"- + ''■<1-'"’ - a'"’) 


(n) 




!_ '• <n> ^*<n>l «P^"^ 

2 . ^ ap<"> ^«i 


-• ^ - .... „*(n) 

+ ■(. 7P^"^) - ^ 


at 




^ ""dp 
dt 


<n) 


+ 7.<L 


(n) 


TP 


<n) 


) - <7 


(n>'.dP 


,*<n).. - 


•' da . 


r I dv 

-i 


+ r TP^"^ - TP*^"^'v f 

J : ? OOL . 

o ^ 


<n> 


t da 


,*<n> 


dt 


' ^ 

da. 

1 


*<n) 


J 


dQ = 


(i = 1,2,.... n> 


0 , 

(6.19) 
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and J " 

V ~ - da . ; t =CD 


dv 


(n) ^<n) _ (n) ^*(n) ” 


+ J TT' "' - L'"' TT 

O 


, -. • dQ = 0 . 

■' da . f t =00 


<i = 1.2 n) , 


( 6 . 20 ) 


In (6.20), the subscript denotes that the parameters , c* . , are 
evaluated at the moment t = oo. Taking into account the 

transversal ity conditions (6.20) and the given initial 

/ 

conditions, we can solve the second order ordinary differential 
equations (6.19) to get the parameters a. and thus the fields 
r and r respectively. 


6.3 FORMULATION OF 6PDP FOR BENARD-DARCT PROBLEM 

The balance equations of mass, momentum and energy for the 
linearized Benard convection in porous media are. 


V»V = 0 . 


dV 
^o dt 


+ V‘P 


ga n T 


t!X 

K 


o C + 7»J 

o V dt q 


C P ft ri 
V o 


•V . 


( 6 . 21 ) 

( 6 . 22 ) 

(6.23) 


where, P is the pressure tensor and is given by , 

o 

= = = VS 

P = p<$ + P 


(6.24) 


=VS 

p being the hydrostatic pressure and P is the symmetrical part 
of the pressure tensor whose trace is zero, n = <0, 0. 1) is 
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the unit vector, V and T are the perturbation velocity and 
temperature fields respectively . denotes the specific heat 

at constant volume and p is the density at temperature T„. J 

^ 0 q 


denotes the 

heat 

current density, g 

i s the 

acceleration due 

to 

gravity, a 

the 

coefficient of volume 

expansion, ft 

the 

temperature 

gradient, t the time 

and 

K represents 

the 

permeabi 1 i ty 

of 

the porous medium. 

<5 i s 

the unit tensor. 

In 


the formulation of Gyarmati’s principle for thermohydrodynami cal 
systems, it is preferable to use the energy picture of the 
principle. In the energy picture, the actual variable is InT 
instead of T and therefore, the balance equations (6.22) and 
<6. 25) take the following form , 

p = p ga n InT - ^ , <6. 25) 
*^0 dt '^o K 


dlnT 

dt 


+ V* 


j = 
q 


Po^ 


<6.26) 


Since we use the energy picture of the principle instead of the 
entropy picture, we replace the entropy production, cr, by the 
energy dissipation. Ter. The energy dissipation for the problem 
is given by , 

o 

=VS o 

To- = - J ‘VlnT - P : <W) . (6.27) 

q 


where 


<W)® 

aft 


1 

2 


dV 

ax 



, <a, /3= 1.2,5) . 


<6.28) 
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The linear laws in this case are 


■^q = ■ ^InT. (6.E9a> 


o 

=VS o 

P = - Lg <W)®. (6.29b) 


where, = X is the coefficient of heat conduction and = Z/u , 
^ being the coefficient of viscosity. The dissipation 
potentials in the energy picture are , 


* 1 ^ 2 s s"^ 

y, = Ty/ = 4- 1_ <71nT) + L (W> : (W> 5 , (6.50) 

2 ^ s J 

o o 

* 1^2 > 

= T(^ = i (J^) + Rg P : P ■ . (6.51) 


Here ^ and ^ . 

The energy picture of Gyarmati's principle is , 


6 J (To- - y>* - 4>*'> dV = 0 . (6.52) 

V 

Principle (6.52), with the help of (6.27), (6.50) and (6.51) 


takes the form , 


=VS o 

6 r - J • VlnT - P : (W) 

J q 

V 

o o , Z 

- \ (VlnT)^ (W)® :(7V)® - 

o o 

, =vs =vs - 

- i_ P ;P r dV = 0 . 

AfJ J 


(6.55) 
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The use of the following vector identiti 


V* (J InT) = InT <7*J ) + j ‘VinT 

M q q • 


=VS 

y. (P .v> 


o o 

=VS =VS o 


y»<p ‘V) = v*(7*p ) + p ;(W)^ 

reduces the principle <6.55) to the form . 


(6.54) 


(6.55) 


<5 J InT 7*J 
V 


^ - J (VlnT)^ - J ^ + v<V*P ) 
H ^ cK q 


O o 1 =VS =VS n 

- fJ (TV)® : (TV)® - 4- P ; P | dV = 0 . 

4 ^/ i 


(6.56) 


The surface integral vanishes due to the boundary conditions on 

o 

=VS 

the two surfaces. Substituting the values of T.J and T.P 

q 

from the balance equations (6.25) and (6.26) and using (6.24), 
we get the principle in the following form , 

6 r (InTXp C /? n*V - p C ) - y (TinT)^ 

J . V o w at c 


- 1 + V* (p g a n InT - p ^ ^ 

2X. q *^0 O ot K 


o o 

O o =VS =VS - 

- p (TV)® :(TV)® “ ^ P 'P ' dV = 0 


(6.57) 


The pressure term vanishes from the volume integral due to the 
boundary conditions. 



V, 

= 0. 

InT = 0 . 


Xj = 0 ; 

5 




Xj = d ; 

^^5 

= 0. 

InT = 0 . 

C^ 

00 
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6. 4 APPLICATION OF DUAL FIELD METHOD 

We introduce a second set of variables In T* and V* which 
are related with the thermodynamic currents through the 
following constitutive equations 


o 

=VS 

P 


71nT . 


- L (W*)® 
s 


<6.59) 


<6. 40) 


Here, the coefficients and are unchanged. It may be 

-M- 'K- 

assumed that InT and V satisfy the same boundary conditions 
as InT and V. These assumed temperature and velocity fields 
are to be determined with the help of balance equations (6.25) 
and <6.26) which become now . 


p - V*<2 aj< 7V)®) + Vp = p ga n InT 

o at o 


- 

K 


<6.41) 




<6.42) 


The principle <6.57) may now be written with the help of <6.59) 
and <6.40) as , 

S J ' <lnT) P n-v - | WlnT)^ 


<VlnT*)^ + V* <p g a n InT - p^ ^ - ^ ) 


9V fjV 




- p <^7V)®:<W)® - p <W ) : <W ) > dV = 0 , 


<6.45) 
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Integrating by parte the terms containing InT* and V* and 
substituting the values of g^nT* and vV 
(6.42). the principle (6.45) reduces to 


then 

and 


6 J <lnT) ip C ft n*V - p C ' - h 

^ o V *^0 V at ’ 


(7inT)' 


1 InT* ip C ft n«V - p C 

2 o V ' at 


+ V* (p^ g a n InT 


P 


o at ~ K 


- p(7v>®:(yv>® 


1 

2 


<P„ g a n InT - p 

fl ■ 


av 

at 


- ^ 
K 


> IdV = 0 


<6.44) 


Here, again the pressure term vanishes due to the boundary 
cond i t i ons . 

To evaluate the principle, we expand each perturbation V 
and In T in complete two-dimensional series of normal modes over 
the plane Xj , x^ , perpendicular to the thermal gradient in the 
following way , 


IX « j ^ ^2^2 at 

InT = ft d 9 <Xj) cos ^ cos ® 


V. 


V «t 

H d d" ® 


V, 


1 dG 

— ^ V . sin 

2 dx- 

a 5 


ajXi 


COS 


2 dG 11 

-4 V j cos . sin 
2 dx, cl 

a ? 


at 

- d ® 


^2^2 wt 
- d ® 


<6.45) 


where, a = <aj + is the wave number of the disturbance 
and a is the frequency which in principle, may be a complex 
number. The velocity components and v^ satisfy the 
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( 


< 


p C V 




'’*4 - 2a^ + sSg* = 3 " e 

dz dz'^ 


+ i- < - a^G + <■ ■'! - a^G . 


Oi: ' . 2 

dz 


^ ^4 2 

dz 


<6.48) 


<6. 49) 


For the linearized Benard convection, the principle of exchange 
of stability is valid, that is, « is real (Chandrasekhar 1961). 
The marginal stability curve can be obtained by setting Real 
<*> = 0 . Since w is real in our case, we put « = 0 in equations 
(6.47), (6.48) and (6.49). Thus the equations of motion and 
GPDP reduce to 


2 

- a )© 


p C V 
o V 


(6.50) 


dz 


dz 


4 


2a^ + a*)G* = ^ ° P d'* a^ ^ 


dz 




, 1 .d" ^2.„ 

+ Dir 2 ^ • 

dz 


(6.51) 


G e - + a^e^i - i G e' 


i. ° 


+ ; p g e. P V G e - Y " g' ^ ^ 

^ ^ .4^ I Z dz J . 

d Da^ a 

Po*" ' , .dG.Z ^ 2^.2 , 1 1 

_4 2 'dF> ^ ° ^ 2> r 

2d t a dz -p 


1 ; 


Pr.'^ 


O i 1 cfG dG 


J 2 p^g a p . G- e - - -GG' 4 ^ gf 5“ p dz = 0 . 

v” dDa'- a •'>- 


(6.52) 
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The equations of motion (6.50> and (6,51), the principle (6.5Z> 
and the boundary conditions <6.47) are quite general for the 
linearized Benard convection in porous media and they are 
applicable to any kind of bounding surfaces. We shall now 
consider the case of two rigid bounding surfaces. 


6.5 SOLUTION FOR RIGID SURFACES 


of rigid 

surfaces , 

the boundary 

conditions are 

0 : G = 

0. 

dz 

0. © = 0 . 


1 ; G = 

0. 

dz 

0. © = 0 . 

<6. 53) 


The similar conditions are to be satisfied by © and G . 
Consistent with the above boundary conditions. we choose very 
simple trial functions , 


2 3 4 

G = A <z - 2z^ + z ) , 


(6.54) 


© = B <z - z^-) . <6. 55) 

where, A and B are the two variational parameters and are to be 
determined from the variational formulation <6.52). Solving the 
equations <6.50) and <6.51) with the help of <6.54) and <6. 55). 
we get © and G as , 


p C V T 4 
• a 


3 2 2 

22 ^ (a'^ + 12)z^ 12z 

2 4 " 4 

a a a 


^ (2.2+24)2^ ^ ^ ^-az 

6 1 2 j 

a 


<6.56) 
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B goi/?d 

2 ' 2 
a 


■^ + 2 - + (v',2 + ¥'x> 


+ <V'^z + e 


+ 2A r . ,2 a^z^ 25 a^z^ 

a Da s- a ^ ^ 


+ (0^2 + 6 


(6.57) 


where . 


2<a^+12> e“® -1 


a -a’ ^2 


2{a^ + 12> 1 - e®; 

6 a -a 

a e -e 


4^5 - ^ 5''2 

b2di - bjd^ 


bjd5 




bid2 - b^dj 


- - 1 - 2 a ¥^4 - ^3 . ¥'6 


^2^4 - ^^2 
b2dV - bjd^ • 


- b,dj 

bjd2 - b^d^ • 


- 6 + a -2a ” ^ 5 • ^6 


“2 + 1 - ^4 . 


/I j. N 2a . , . 2a . , 2 

(1 + a) e + a - 1 , b^ = ae + 2a - a . 


e - a + 5 , 


b, = 6e® + a^ - 6a + 18 
4 


j 2a , 

dj = e - 1 , 


d^ = e"” - 2a - 1 . 


d, =— ^<e -1)+1-— . 

5 a^ ^ 


. 12. a... a 12. _ 

d. =— y<e-l>+e -a-— +5. 
4 2 a 

a 
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The expressions <6.56> and (6.57> satisfy the 
cond i t i ons . 


z = 

0: G = 0, 

55*- 0 

dz 

©* = 0 , 


z = 

1 : G ~ 0 , 

0 

dz 

9 * = 0 . 


Final ly 

the GPDP 

(6.52) , 

with the 

subst i tut i on 


expressions for G. 9 * and G* . becomes . 


6 r p c V _ AB _ 1 2 

i^ov ' 140 2 




r I., 2 AO P ^ A 1 o 

^ '’o ^4 140 -^4 63b - ^2 


'’o- 4 ^ a' ^ 4 , ^2 

2^4 105 630 5^2 > * 


kv‘ 


Rk B 


2 AB 


A R < ^ ^ /?, > 


2 *^0 .4 

d 


2 ' 2 ■ 4-, ''3 

a a Da 


+ 1 < Rk AB « + 2A^ « > 

2 .4_,^ ^ V 2 ^4 ^ 4^ ^5 ^ 
d Da a a Da 


^ i '=’0 ^ 2Rk ^ 4A^ 

2 .4 „ 2 ^ V 2 ^6 ^ 4r, ^7 ^ 

d Da a a a Da 


0 , 


where, R = " ^Rayleigh number), k = 

^o'^v 


a" - 3a2 




105 630 ; 


i + ¥', <C, - 2C, + C„) 


- ^ ‘'10 


boundary 

(6.58) 

of the 


(6.59) 


* '"2 “=6 
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^4 


ft 


6 


^7 


50 " 3 ^ + + ^4 




+ V', 

5 ^^6 

- 

+ 

. <<=4 - 


1 

2 

2 

a 






■ To 

2 

a 

280 

^3 <c 

2 ■ ^7 

) + 

*4 

- S' 




+ <P, 

- 

s> 

♦ *6 <'=4 - ■ 

_ 1 

2 


<s - 





140 

15 a 

2 " V ^3 

2C5 + 

^11 

> + ¥^4 

<C^ - 2C^+ C^) 


+ V '3 

- 

^^10 ■ 

^ ^12> 

+ ^ 

6 <^6 - 

+ So' • 

_ 1 
105 

2 

5 a^ 

2 

a 

1260 

+ ^3 

<S - 

'S 

" ‘^ll^ 


+ ^4 

<s ■ 

- 2C7 + 

S> + ^5 (Cj - 

2^10 + 

S2' 


" ^6 

- 

2C8 + 

^10^ ’ 





■ 50 ■*■ ^^^3 " 5 Cj + 20^) + a (Cj - 30 ^ + 2 C^>} 

+ a <Cj - 3 C^ + 2C7) 

+ ¥-5 C<c^ - 3C^ + 2Cg> - a (C^ - JCg + 2Cj(,>} 

- a <c^ - 5 C^ + 2c^, , 

1 _ a^ 

5 210 ■•■ *■} ^‘■'5 ■ 5 Cj + 2C^) + a <Cj - 5C^ t ZC,)! 

+ a (Cj - 3 Cj + 2C7) 

+ <>5 C<C^ - 3 C^ + ZCj) - a (Cj - 5 Cj + 2 Cj (,>3 

- ’Cj + 2Cj> , 


- a v-j (C^ 



14 ] 


1 


0 


1 


= s dz = 


1-e 


-a 


a 


= J dz = i! - 

a a • 

0 

1 

= f ze-^^ dz = - t fi 

a a ' 

0 


.2_az . ^^5 


r az _■ e 
- J z e dz = I- - 


0 

1 


a a 


2.-az , e“® ^ ^^4 


a a 


'10 


= / dz = 

0 

= / dz = . 

a a ' 

0 

1 

r 5 ~az e~^ 

= \ z e dz = - - + - ^ 

a a 

0 

= jVe- dz = l! - , 

•' a a 

0 

= /z*e-dz 

^ a a 


- r H ^*"9 

-jl - J z e dz = j ^ , 

0 

hz -s = - C " 
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Taking variation of (6.59), the marginal stability curve is 
obtained as , 

a^<a^+10> (1 + 1^) 

« 

+ ) . 

numbers and the corresponding 
are obtained for various values of 
in Table 1 . 

Table 1 Critical values of a and DaR 
Da a DaR 



00 (Nav i er-Stokes ) 5.12 


1748.697 ( 

= R ) 
c 


lO"^ 

5.50 


44.794 



10"^ 

5.51 


42.557 



lo"^ 

5.51 


42.555 


The 

results 

given in Table 1 show 

that 

when the permeabi 

lity of 

the 

medium 

tends to infinity. 

the 

system behaves 

like an 

ord i 

nary fluid, with critical wave and Rayle igh-Darcy 

number 


equal to 5.12 and 1748.697 respectively. These two values are 
remarkably close to the exact values of 5.117 and 1707.762 given 
by Chandrasekhar (1961). Results obtained for different values 


R 


ua 


'50" 


140 ^4®- 

where. E = (-/3^ + \ 

„4 ' 5 2 

a 


6.6 RESULTS AND DISCUSSION 


The critical wave 
Ray 1 e i gh-Darcy numbers (DaR) 
Da. These results are shown 
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of Ray 1 e igh-Darcy number are also comparable with the critical 
wave number and corresponding Ray 1 e igh-Darcy number determined 
by Lebon and Cloot <1986) which differ by less than 3* . This 
shows that the governing principle of dissipative processes 
proposed by Gyarmati is well suited for analyzing the Benard 
convection in porous media. 
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